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PREFACE 


This  report  is  the  first  in  a series  dealing  with  ground  control  for 
tunnels  constructed  by  shield  techniques  in  soil.  It  addresses  the  issue  of 
the  contribution  of  time-dependent  dissipation  of  excess  pore  pressures  to 
the  long-term  behavior  of  the  ground  around  a tunnel  and  the  liner  support 
system.  Because  most  tunnel  monitoring  to  date  has  concerned  only  short- 
term response,  this  work  uses  primarily  analytical  methods  to  generate  the 
necessary  data.  However,  the  available  field  information  is  also  reviewed. 

Other  reports  to  be  issued  in  the  series  associated  with  this  one 
focus  exclusively  on  ground  control  obtained  by  using  advanced  shield 
tunneling,  a subject  also  briefly  addressed  herein.  The  research  sponsor- 
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for  the  project. 
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EXECUTIVE  SUMMARY 


A major  technical  challenge  with  shallow  tunneling  in  urban,' areas 
is  to  limit  settlements  of  the  ground  surface  which  may  cause  damage  to 
buildings.  The  question  of  how  much  settlement  will  occur  has  largely 
been  studied  empirically  in  the  past.  However,  as  refinements  in 
settlement  control  and  construction  techniques  are  developed  it  becomes 
more  important  to  understand  the  mechanisms  of  these  problems. 

One  particular  aspect  of  ground  response  to  tunneling  which  is  only 
poorly  understood  is  the  effect  of  mobilization  and  dissipation  of 
excess  pore  pressures  in  cohesive  soil.  Based  on  a limited  amount  of 
field  data,  it  would  appear  that  the  delayed  settlements  which  are 
related  to  pore  pressure  dissipation  can  be  larger  than  those  which 
occur  immediately  after  the  shield  passes.  The  available  data  suggests 
that  time-dependent  settlements  can  be  a problem  for  both  conventional 
as  well  as  advanced  shields.  Time-dependent  movements  are  also  noted  to 
be  strongest  where  the  shield  procedures  lead  to  formation  of  a large 
tail  void  or  otherwise  cause  significant  remolding  of  the  soil.  In  the 
case  of  advanced  shields  using  a pressurized  face,  consolidation  effects 
may  be  influenced  by  initial  outward  heaving  of  the  soil  at  the  face. 

In  addition  to  time-dependent  settlements,  a number  of  field 
measurements  have  shown  that  long  term  changes  in  liner  diameter  and 
load  distribution  on  the  liner  occur.  The  observations  suggest  that  the 


radial  load  and  horizontal 


liner  diameter  increase  approximate! y 


linearally  with  the  logarithm  of  elapsed  time  after  construction.  This 
behavior  has  important  implications  for  liner  design  as  well  for  repair 
on  construction  alteration  work  undertaken  some  years  after  the  original 
tunnel  construction  is  completed. 

This  investigation  was  undertaken  to  study  those  aspects  of  the 
time  dependent  tunnel  response  which  are  due  to  dissipation  of  excess 
pore  pressures  generated  in  the  ground  during  tunneling.  For  this 
purpose,  a finite  element  code  is  developed  which  simulates  tunneling  in 
soil,  explicitly  taking  into  account  pore  pressure  mobilization  and 
dissipation  in  time.  The  program  provides  a means  of  investigating  time 
dependent  effects  on  tunneling  settlements,  and  liner  loads  as  related 
to  consolidation  in  the  soil  around  the  tunnel. 

The  program  incorporates  an  e 1 asto-p 1 ast i c soil  model  and  large 
strain  effects  in  order  to  reasonably  simulate  the  behavior  of  the  soil 
as  it  flows  plastically  to  fill  the  tail  void  around  the  tunnel.  Other 
aspects  of  the  tunneling  problem  are  simulated  as  realistically  as 
possible.  Provisions  are  made  for  modelling  gravity  stress  profiles, 
soil  properties  varying  with  depth,  and  liner  installation  after 
mobilization  of  a user  specified  degree  of  soil  closure  into  the  tunnel. 

Parametric  studies  on  tunnel  construction  in  clay  are  performed 
into  the  effects  of  soil  strength,  tunnel  liner  stiffness,  size  of  the 
tail  void,  initial  soil  heave  and  permeability  of  the  soil.  The  results 
indicate  that  time-related  consolidation  effects  around  tunnels  in 
cohesive  soil  can  be  of  great  importance.  Surface  settlements  and  liner 
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loads  generally  increase  with  time  until  the  consolidation  process  is 


completed.  The  increase  in  settlements  is  observed  to  be  large  when 
flexible  liners  are  used,  soft  soils  are  present,  or  the  tail  void  is 
large.  Small  degrees  of  initial  outward  heave  due  to  a pressurized  face 
shield  seem  to  have  beneficial  effects  since  settlements  are  reduced 
relative  to  cases  where  only  inward  movement  occurs. 

The  result  of  the  parametric  studies  also  show  that  liner  loads, 
load  distribution  and  diameter  change  with  time.  Typically  the 
predicted  trends  are  similar  to  those  observed  in  the  field  during  the 
consolidation  process.  With  time,  the  load  distribution  tends  to 
accentuate  the  bending  stresses  in  the  liner  and  the  total  loads 
increase  slightly.  The  increases  in  stresses  in  the  liner  vary 
approximately  linearally  with  the  logarithm  of  time  until  the  excess 
pore  pressures  are  dissipated.  A similar  variation  is  observed  for  the 
lengthening  of  the  horizontal  diameter  with  time.  These  results  in  part 
explain  the  heretofore  little  understood  changes  which  have  been 
observed  in  the  field. 


Chapter  1 


INTRODUCTION 

Underground  construction  in  urban  areas  is  needed  for  rapid  transit, 
water  supply  and  sewage  disposal.  In  the  past,  underground  construction 
in  the  United  States  was  usually  performed  using  cut  and  cover 
techniques.  Today,  however,  tunneling  is  becoming  more  popular  because 
of  recent  advances  in  shield  tunneling  technology  and  because  it  causes 
less  environmental  disruption. 

One  major  technical  challenge  with  shallow  tunneling  in  urban  areas 
is  to  limit  settlements  of  the  ground  surface  which  may  cause  damage  to 
buildings.  The  question  of  how  much  settlement  will  occur  has  largely 
been  studied  empirically  in  the  past.  However,  as  refinements  in 
settlement  control  and  construction  techniques  are  developed  it  becomes 
more  important  to  understand  the  mechanisms  of  these  problems. 

One  particular  aspect  of  ground  response  to  tunneling  and  related 
activities  which  is  only  poorly  understood  is  the  effect  of  mobilization 
and  dissipation  of  excess  pore  pressures  in  cohesive  soils.  In  his 
excellent  state-of-the-art  paper,  Peck  (1969)  noted,  "Delayed 
settlements  due  to  long-time  consolidation  of  the  clay  around  the  tunnel 
should  not  be  overlooked.  These  may  spread  much  more  widely  than  those 
due  to  the  tunneling  operations  themselves."  This  behavior  has  become 
even  more  significant  where  advanced  shield  techniques  are  used,  since 
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there  is  field  evidence  that  the  soil  in  this  case  is  heaved  aside  as 


the  shield  moves  forward  (Clough,  et  al.,  1982).  The  soil  around  the 
shield  is  thus  remolded,  and  time  dependent  movements  are  observed  due 
to  consolidation. 

Probably  the  main  reason  for  a lack  of  study  of  the  consolidation 
problem  is  that  it  is  complex  and  difficult  to  characterize 
analytically.  Also  few  field  data  are  available,  since  most  actual 
field  measurements  are  made  only  during  or  immediately  after 
construction,  and  pore  pressures,  per  se,  are  rarely  monitored  at  all. 

This  report  presents  to  the  knowledge  of  the  authors  the  first 
comprehensive  investigation  of  time-dependent  effects  due  to  the 
dissipation  of  excess  pore  pressures  around  tunnels.  In  particular,  the 
following  items  are  addressed: 

1.  Available  field  documentation  on  time  effects. 

2.  Development  of  a finite  element  code  which  allows  simulation  of 
tunneling  in  cohesive  soil,  explicitly  taking  into  account  pore 
pressure  mobilization  and  dissipation  in  time. 

3.  The  mechanism  of  excess  pore  pressure  development  and  dissipation 
around  tunnels  in  cohesive  soil  and  its  effects  on  time-dependent 
settlements  and  liner  loads. 

4.  The  mechanism  of  excess  pore  pressure  development  and  dissipation 
in  cohesive  soil  during  advanced  shield  tunneling  with  a 
pressurized  face. 
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From  the  results  of  these  studies  an  assessment  is  developed  of  the 
significance  of  time-related  consolidation  effects  around  tunnels  in 
cohesive  soil. 

In  order  to  set  the  stage  for  the  analytical  studies.  Chapter  II 
presents  a review  of  the  existing  literature  on  time-dependent  tunnel 
behavior.  Also,  a technique  for  realistically  simulating  the 
construction  procedure  of  a tunnel  using  the  finite  element  method  is 
introduced.  This  includes  modelling  in  plane  strain  the  excavation  of 
the  soil,  the  installation  of  a liner  and  the  closing  of  a non-uniform 
tail  void  gap. 

The  finite  element  method  is  used  because  it  allows  the  coupled 
equations  of  elasticity  and  fluid  flow  to  be  solved  for  complicated 
boundary  and  soil  conditions.  Closed  form  analytical  solutions  to  such 
problems  are  impossible.  The  finite  element  method  is  particularly 
useful  for  performing  parametric  studies  which  may  be  used  to  identify 
the  most  important  factors  contributing  to  the  time-dependent  behavior 
of  surface  settlements  and  liner  loads. 

An  e 1 asto-p 1 ast i c soil  model  is  needed  to  provide  a good 
representation  of  the  stress-strain  behavior  of  the  clay,  since  it 
undergoes  yielding  as  it  moves  into  the  tail  void.  The  Cam  clay  soil 
model,  an  e 1 asto-p 1 ast i c strain  hardening  model  is  used  herein.  This 
model  was  developed  at  Cambridge  University,  where  it  has  been  used  in 
both  experimental  and  analytical  tunnel  modelling.  A method  for 
implementing  it  in  a two-  or  three-dimensional  code  is  presented  in 
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Chapter  III. 


A theory  for  e 1 as to-p 1 as t i c consolidation  is  outlined  in  Chapter  IV. 


This  theory  allows  the  equations  governing  the  ground  movement  to  be 
coupled  with  those  governing  the  flow  of  water  through  the  soil  voids  in 
order  to  calculate  the  excess  pore  pressures.  A pore  pressure  degree  of 
freedom  is  added  to  some  of  the  nodal  points.  In  order  to  incorporate 
this  theory  into  the  finite  element  code,  a numerical  stability  study 
was  used  to  choose  a suitable  finite  element  and  integration  rule  for 
the  time  domain.  Using  this  method,  the  development  and  dissipation  of 
pore  pressures  in  time  may  be  modelled  for  the  case  of  a tunnel  in  a 
saturated  cohesive  soil  typical  of  many  urban  areas. 

Large  strains  occur  near  the  tail  void  as  it  is  closed.  Since  this 
closing  is  the  major  cause  of  the  surface  settlements,  it  is  important 
that  a large  strain  formulation  be  used.  A simple  approach  is 
demonstrated  in  Chapter  V,  which  gives  accurate  solutions  even  when 
displacements  are  large.  The  method  is  efficient  since  the  symmetry  of 
the  stiffness  matrix  is  maintained. 

The  final  finite  element  code  which  incorporates  all  of  the  options 
offers  a tool  which  will  accurately  simulate  the  construction  of  a 
tunnel  including  predictions  of  surface  settlements,  liner  stresses  and 
excess  pore  pressures.  Results  of  parametric  studies  on  tunnel 
construction  in  clay  are  given  in  Chapter  VI.  The  effects  of  varying 
soil  strength,  tunnel  liner  stiffness,  size  of  the  shield  tail  void,  and 
permeability  of  the  soil  are  considered.  Also,  the  influence  of  an 
initial  outward  movement  of  the  soil  is  examined.  This  effect  is 
modelled  since  heaving  can  occur  during  advanced  shield  tunneling. 
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Finally, 


comparisons  are  made  between  finite  element  behavior  and 


existing  field  data. 

The  final  chapter  provides  a summary  of  the  work  done,  and 
conclusions  and  recommendations  are  given.  Following  this  chapter  three 
appendices  are  included  which  describe  the  computer  code  and  present  a 
User's  Guide  for  it. 
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Chapter  2 


TIME-DEPENDENT  BEHAVIOR  OF  SHIELD  DRIVEN  TUNNELS 
2 . 1 INTRODUCTION 

In  this  chapter  both  existing  field  data  and  theoretical  studies  on 
the  time-dependent  behavior  of  tunnels  are  reviewed.  A method  for 
modelling  the  construction  procedure  and  subsequent  time-dependent 
behavior  of  shield  driven  tunnels  using  finite  elements  is  given. 
First,  current  technology  in  shield  tunneling  methods  is,outlined. 

2.2  REVIEW  OF  SHIELD  TUNNELING 

Tunnels  in  soft  clay  are  usually  constructed  with  the  aid  of  a 
shield.  The  shield  is  basically  a steel  cylinder  which  provides 
temporary  support  for  the  soil  while  excavation  proceeds  at  the  tunnel 
face,  as  shown  in  Figure  2.1.  Liner  segments  are  erected  in  the  tail  of 
the  shield.  The  tunnel  is  advanced  by  jacking  the  shield  ahead  while 
thrusting  against  the  in-place  liner. 

A tail  void  is  created  between  the  liner  and  the  soil.  The  tail 
void  size  is  determined  by  the  thickness  of  the  shield  plus  the 
tolerance  gap  between  the  outside  of  the  liner  and  the  inside  of  the 
shield.  The  tail  void  is  further  increased  when  overcutters  are  used  to 
aid  steering.  Plowing  occurs  when  proper  alignment  with  grade  is  not 
maintained  and  this  can  add  dramatically  to  the  size  of  the  tail  void, 
as  shown  in  Figure  2.2. 
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Attempts  are  usually  made  to  fill  the  tail  void  by  pumping  some  form 
of  grout  or  pea  gravel  through  holes  in  the  liner  segments.  In  weak 
soils,  however,  this  void  usually  closes  before  any  material  can  be 
injected.  When  grout  is  being  injected,  the  gap  between  the  shield  and 
the  liner  must  be  large  enough  to  house  grout  seals.  This  may  also 
increase  the  size  of  the  tail  void. 

In  weak  cohesive  soil  the  face  of  the  tunnel  usually  must  be 
supported  during  the  excavation  procedure.  Many  different  techniques 
for  support  exist,  and  in  poor  soil  combinations  of  these  may  be 
required.  Richardson  and  Mayo  (1975)  provide  a useful  review  of 
conventional  procedures  for  this  purpose.  The  most  prominent  technique 
in  the  U.S.  in  the  past  involves  the  application  of  compressed  air  in 
the  tunnel.  Insurance  costs  have  risen  dramatically  when  compressed  air 
is  used  however,  and  alternative  techniques,  particularly  advanced 
closed  face  shields  are  becoming  more  attractive.  Predecesors  to  the 
advanced  shields  were  developed  for  use  on  the  Bart  system  (Peterson  and 
Frobenius,  1971);  this  type  of  shield  had  a closed  face  with  a rotating 
cutter  wheel.  The  cutter  wheel  had  provisions  for  windows  or  doors 
which  were  opened  to  bring  the  soil  in.  However,  where  zones  of  running 
or  flowing  ground  are  encountered  within  the  cohesive  soil,  compressed 
air  has  to  be  used  to  keep  soil  from  moving  freely  through  the  cutter 
wheel  openings.  This  essentially  negates  the  advantages  of  such  a 
mach i ne . 

To  overcome  the  problems  with  compressed  air,  a new  technology  has 
been  introduced  where  the  tunnel  face  is  supported  either  by  a slurry 
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under  pressure  or  by  a combination  of  earth  and  water  pressure.  Clough 


(1980,  1981)  has  provided  a general  description  of  these  shields.  A 
typical  earth  pressure  balance  shield,  is  shown  in  Figure  2.3.  The 
first  application  of  this  shield  in  the  U.S.  took  place  in  1981  in  San 
Francisco  for  a sewer  tunnel.  The  performance  for  this  project  is 
described  in  a companion  volume  for  this  research  project  (Clough,  et 
al.,  1982). 

In  spite  of  the  innovations  in  shield  technology  for  control  of  the 
soil  at  the  face  of  the  tunnel,  a key  problem  remains:  the  tail  void  and 
the  soil  closure  into  it.  Also  excessive  initial  heaving  and  subsequent 
settlement  can  occur  if  the  shield  is  not  operated  properly. 

Surface  settlements  for  shield  tunneling  may  be  divided  into  four 
stages,  (Cording  and  Hansmire,  1975).  They  develop 

1.  ahead  of  the  face. 

2.  over  the  shield. 

3.  during  the  erection  of  the  lining  at  the  tail  of  the  shield. 

4.  with  time  and  further  advance  of  the  tunnel,  as  the  lining 
deflects. 

The  first  three  stages  of  surface  settlement  are  shown  schematically  in 
Figure  2.1.  If  an  advanced  shield  is  being  used,  settlements  ahead  of 
the  face  should  be  small  since  the  in-situ  stresses  at  the  face  will  be 
balanced.  The  settlement  over  the  shield  is  somewhat  larger,  but  the 
settlement  is  most  noticeable  on  the  surface  as  the  tail  of  the  shield 
passes . 
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Figure  2.3:  Typical  Earth  Pressure  Balance  Shield 


Thus,  surface  settlements  are  highly  dependent  on  the  details  of  the 


shield  design  and  the  experience  of  the  operators. 

2 . 3 CAUSES  OF  TIME-DEPENDENT  BEHAVIOR 

Tunnel  behavior  in  clay  has  been  found  to  be  time-dependent.  Both 
surface  settlements  and  liner  stresses  generally  increase  with  time. 
There  are  four  major  causes  of  this  time-dependent  behavior. 

First,  a tunnel  is  constructed  by  advancing  a heading.  A 
three-dimensional  state  of  stress  exists  near  the  heading,  while  away 
from  the  heading  the  state  of  stress  may  be  assumed  to  be  plane  strain 
with  zero  strain  in  a direction  parallel  to  the  axis  of  the  tunnel.  As 
the  heading  advances  with  time  away  from  a fixed  point,  the 
three-dimensional  effect  at  this  point  diminishes. 

Finite  element  techniques  have  been  used  to  investigate  this  type  of 
time-dependent  behavior.  Ghaboussi  and  Gioda  (1977)  studied  the 
time-dependent  effect  of  advancing  tunnels  using  axisymmetric 
conditions.  This  assumption  leads  to  an  approximate  solution  which  may 
give  useful  results  for  deep  tunnels.  Sakurai  (1978)  also  investigated 
this  problem,  reducing  it  to  plane  strain  by  introducing  the  'equivalent 
initial  stress'.  While  this  method  appears  to  give  some  useful  results, 
it  remains  an  approximation  to  the  true  three  dimensional  behavior.  The 
first  full  three-dimensional  analysis  of  advancing  shield  tunnels  is 
described  in  a report  for  this  research  project  (Kasai i and  Clough, 
1982) . 
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A second  cause  of  time-dependent  behavior  is  due  to  the  effect  of 


soil  creeping  or  undergoing  secondary  compression.  Finite  element 
techniques  have  been  used  to  investigate  this  phenomenon  using 
visco-elastic  models  by  Ghaboussi  and  Gioda  (1977),  Sakurai  (1978),  and 
Tan  and  Clough  (1980).  This  effect  was  found  to  be  of  particular 
importance  where  chemical  grouting  techniques  are  used. 

For  tunnels  located  below  the  water  table  there  is  a third  cause  of 
time-dependent  behavior.  In  the  long  term,  the  tunnel  may  act  as  a pore 
water  sink.  As  the  water  drains  into  the  tunnel,  the  water  table  may  be 
lowered,  producing  further  settlements.  Tunnel  liners  are  usually 
water-proofed  in  order  to  prevent  this  effect.  Nevertheless,  it  is 
difficult  to  make  tunnel  liners  totally  impervious.  Fitzpatrick  et  al . 
(1981)  used  the  finite  element  method  to  investigate  the  effect  of  liner 
leakage.  An  elastic  soil  model  was  used  and  the  construction  process 
itself  was  not  modelled.  The  final  settlement  due  to  liner  leakage  was 
found  to  be  proportional  to  the  depth  of  the  tunnel  and  the  tunnel 
diameter.  They  also  found  that  the  flow  of  water  into  the  tunnel  was 
negligible  when  the  liner  permeability  was  two  orders  of  magnitude 
greater  than  that  of  the  soil. 

The  fourth  cause  of  time-dependent  behavior  is  due  to  consolidation 
of  the  remolded  layer  around  the  tunnel,  the  behavior  which  is  the 
subject  of  this  report.  During  the  construction  procedure  there  is  a 
redistribution  of  stress  in  the  neighborhood  of  the  tunnel.  This  stress 
redistribution  sets  up  excess  pore  pressures  in  the  clay.  Further 
settlements  and  changes  in  liner  loads  take  place  as  these  excess  pore 
pressures  dissipate. 
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2 . 4 OBSERVED  TIME  EFFECTS  DUE  TO  PORE  PRESSURE  DISSIPATION 


2.4.1  Spttl ements 

One  of  the  first  documented  case  histories  of  shield  tunneling  in 
clay  produced  evidence  of  significant  consolidation  effects.  The  tunnel 
was  constructed  in  Chicago  in  the  early  1940's,  and  the  performance  has 
been  described  by  Terzaghi  (1942).  The  shield  technique  used  was 
somewhat  unusual  in  that  it  involved  shoving  a closed  face  shield 
through  the  clay  with  small  doors  in  the  face  which  were  opened  to  allow 
the  clay  to  squeeze  through  during  the  advance.  The  clay  was  generally 
soft  to  medium  along  the  tunnel  alignment  with  compressive  strengths 
ranging  from  600  to  1400  psf,  (28.8  to  67.1  kN/m2). 

It  was  observed  that  as  the  shield  moved  through  the  soil,  a 
pronounced  heave  surface  occurred  initially  with  movements,  of  one  to 
four  inches  (2.5-10  cm).  This  apparently  was  caused  by  the  fact  that 
the  volume  of  the  soil  entering  the  shield  through  the  face  doors  was 
less  than  the  advance  volume  occupied  by  the  shield.  Subsequent  to  the 
passage  of  the  shield,  and  for  many  days  thereafter,  the  ground  surface 
settled.  Figure  2.4  depicts  the  time  dependent  nature  of  the  surface 
settlements.  At  this  particular  location  time-dependent  settlements 
were  observed  up  to  a period  of  150  days,  reaching  values  of  eight  to  10 
inches  (20-25.4  cm).  Terzaghi  (1942)  attributed  the  initial  subsidence 
effects  after  shield  passage  to  a partial  collapse  of  the  clay  onto  the 
liner  plate.  But  the  significant  subsequent  time-dependent  movements  he 
felt  were  due  to  consolidation  effects. 
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Figure  2.4:  Observed  Settlement  Versus  Time  for  Tunnels  in  Clay  and 

Silts  in  Chicago  and  Osaka 
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Terzaghi  (1942)  explained  the  consolidation  phenomenon  as  occurring 


because  of  excess  pore  pressures  set  up  in  a remolded  zone  in  the  clay 
around  the  tunnel,  created  as  a result  of  the  heaving  of  the  soil  during 
a shove.  He  noted  that  the  pore  pressures  likely  drained  into  the  sand 
layer  which  was  backfilled  behind  the  liner  plate  after  the  liner 
segments  emerged  from  the  tail  of  the  shield.  It  is  significant  that 
the  settlements  which  ultimately  occurred  were  always  greater  than  the 
initial  observed  heave  values. 

Because  of  the  significant  ground  movements  induced  by  the  shoving 
process  discussed  by  Terzaghi  (1942),  this  technique  has  been  largely 
avoided  in  recent  times.  Open-face  shields  were  preferred  and  largely 
used  until  the  recent  advent  of  the  rotating  cutter  head  machines  with 
pressurized  faces.  In  soft  soils,  compressed  air  has  been 
conventionally  used  with  open-face  shields  to  control  free  stability. 
Few  instances  of  time  dependent  movements  have  been  reported  for 
conventional  or  advanced  shields,  until  the  last  few  years. 

Time-dependent  settlement  above  a sewer  tunnel  in  Thunder  Bay, 
Ontario  has  been  described  by  Belshaw  and  Palmer  (1978)  and  Palmer  and 
Belshaw  (1980).  A rotating  cutterhead  machine  was  used  with  a segmented 
concrete  liner  erected  within  the  tail  of  the  shield.  The  excavated 
diameter  of  the  shield  was  8.1  feet  (2.47  m)  at  a depth  of  31  feet  (9.5 
m)  to  the  crown.  The  tunnel  was  constructed  in  a soft  clay  with  a shear 
strengh  of  730  psf  (35  kN/m2)  at  the  tunnel  springline.  This  clay 
extended  from  a depth  of  about  20  feet  (6  m)  to  the  bedrock  at  a depth 
of  80  feet  (25  m).  A layer  of  silt  lay  above  the  soft  clay  and  the 
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water  table  was  located  at  a depth  of  5 feet  (1.5  m). 


The  surface 


settlement  above  the  tunnel  centerline  was  monitored  at  a test  section, 
and  was  found  to  increase  rapidly  as  the  tail  of  the  shield  passed  the 
test  section.  After  about  12  days  the  settlement  appeared  to  stop 
increasing.  The  settlement  then  began  to  increase  again  and  was 
monitored  for  one  year.  This  settlement  is  plotted  versus  time  in 
Figure  2.5.  After  one  year  the  settlement  had  increased  by  100  percent, 
a behavior  attributed  to  consolidation  effects  as  the  excess  pore 
pressures  in  the  region  of  remolded  clay  dissipated.  Little  to  no 
drainage  occurred  into  the  tunnel  since  a clay  seal  was  grouted  in 
around  it  and  only  a minor  drawdown  was  noted  in  piezometers  located 
near  the  tunnel.  It  may  be  noted  also  that  there  is  a strong  similarity 
between  the  form  of  the  results  in  Figure  2.5  and  those  reported  by 
Terzaghi  (19-12)  for  the  tunnel  in  Chicago. 

The  behavior  of  a later  section  of  the  tunnel  in  Thunder  Bay  was  also 
monitored.  Here  the  settlement  was  found  to  increase  by  only  15  percent 
during  a similar  time  period  as  observed  at  the  first  test  section.  The 
difference  in  behavior  was  attributed  to  the  greater  clay  cover  in  the 
earlier  section  and  to  shield  alignment  difficulties  during  the  early 
driving  stages  at  the  first  test  section.  Where  the  problems  with 
shield  steering  occurred,  presumably  a greater  zone  of  remolded  clay  was 
produced,  giving  rise  to  large  future  settlements.  This  is  an  important 
finding  and  suggests  that  quality  of  construction  procedures  has  a large 
effect  on  the  amount  of  consolidation  after  shield  passage. 
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Figure  2.5:  Observed  Settlement  Versus  Time  for  a Tunnel  in  Clay 


The  most  recent  incidents  of  time  dependent  ground  movements  reported 


in  the  literature  have  been  for  earth  balance  shield  projects. 
Kitamura,  et  al.,  (1981)  described  the  observed  surface  movements,  which 
occurred  for  a 22.1  ft.  (6.7  m)  diameter  shield  tunneling  in  silts  and 
clays  at  a depth  of  15-22  ft.,  (4.5-12.7  m).  The  movements  are  shown  in 
Figure  2.4;  as  the  shield  approached  the  instrumentation  station,  a 
heave  of  0.5  in  (1.3  cm)  was  measured.  As  the  shield  passed  the 
instrumentation  line,  settlements  began  to  occur,  and  they  continued  for 
a period  of  150  days  (see  Figure  2.4). 

The  initial  rapid  settlements  immediately  after  shield  passage  are 
likely  due  to  the  tail  void  effect.  However,  there  is  an  extended 
delayed  settlement  which  eventually  leads  to  a doubling  of  the  initial 
movements.  Importantly,  the  form  of  the  settlement  pattern  is  very 
similar  to  that  described  previously  for  the  projects  in  Chicago  and 
Thunder  Bay.  The  delayed  movements  would  appear  to  be  related  to 
consolidation  of  the  soils  around  the  tunnel.  Presumably  excess  pore 
pressures  were  set  up  by  the  heaving  effect  in  front  of  the  shield  and 
the  inward  movement  toward  the  tail  void. 

The  second  earth  pressure  balance  shield  project  where  long-term 
settlements  occurred  was  instrumented  and  monitored  in  1981  by  a team  at 
Stanford  University  as  a companion  effort  to  the  work  reported  herein. 
Details  of  the  results  are  given  elsewhere  by  Clough  et  al . (1982),  and 
only  the  time-dependent  response  is  described  for  this  work.  The 
tunnel,  12  ft.  (3.6  m)  in  diameter  and  at  a depth  of  about  30  ft.  (9.1 
m),  is  located  in  San  Francisco.  Along  the  alignment  the  soil  profile 
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consists  of  about  20  ft.  (6.1  m)  of  rubble  fill  underlain  by  soft  silts 


and  clays.  The  uatertable  is  at  a depth  of  about  10  ft  (3  m). 

Both  vertical  and  lateral  movements  were  measured  during  and 
following  shield  passage  at  three  instrumentation  stations.  The  data 
shoued  that  most  of  the  movements  were  concentrated  in  a zone  around  the 
shield  and  in  primarily  the  soft  soils.  The  rubble  fill  apparently 
acted  as  a structural  element  and  smoothed  the  deflections  over  a 
broader  area.  The  lateral  movement  devices  shoued  very  clearly  that  the 
shield  initially  heaved  the  soil  aside  and  away  from  the  face.  While 
very  little  initial  movement  was  observed  at  the  surface,  up  to  three 
inches  (7.6  cm)  of  outward  displacement  occurred  in  the  soft  soil  around 
the  shield  during  each  shove.  Following  shield  passage,  the  soil  began 
to  move  inwards  towards  the  tunnel,  and  this  trend  has  continued  to  the 
time  of  this  writing  (December  1981).  Figure  2.6  shows  a plot  of 
measured  surface  movements  versus  time  at  one  instrumentation  station. 
As  was  the  case  with  the  earth  balance  shield  behavior  reported  by 
Kitamura,  et  al . (1981),  there  is  an  initial  heave  during  shield 
passage.  This  is  followed  by  settlement  which  occurs  over  a period  of 
150  days.  Small  levels  of  continuing  movements  are  continuing  to  be 
observed  in  the  instrumentation  at  this  site  as  of  this  date  (December, 
1981)  and  it  is  intended  to  continuing  the  monitoring  process  for  at 
least  one  year. 

In  summary,  the  review  of  observed  ground  movements  has  shown  that 
significant  time  dependent  effects  are  common.  This  is  true  for  both 
conventional  and  advanced  shields. 
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Figure  2.6:  Observed  Ground  Movements,  Versus  Time  For  An  Earth 

Balance  Shield  in  San  Francisco  (After  Clough,  et  al . , 
1982) 
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The  degree  of  the  time  dependent  movements  also  appears  to  be  a function 
of  the  degree  of  disturbance  produced  by  the  tunneling  techniques.  The 
more  the  disturbance,  the  larger  the  time-dependent  movements.  All  of 
the  information  available  suggests  that  the  movements  are,  in  large 
part,  due  to  consolidation  phenomena  related  to  the  dissipation  of 
excess  pore  pressures. 

2.4.2  Liner  Loads  and  Shapes 

The  load  acting  on  a liner  and  the  horizontal  diameter  of  the  liner 
have  been  observed  to  generally  increase  with  time.  Peck  (1969).  Many 
field  measurements  of  this  type  have  been  reported  for  tunnels  in  London 
clay,  and  limited  results  have  been  reported  for  other  soils. 

Figures  2.7  and  2.8  show  typical  increases  in  liner  loads  and 

© 

horizontal  tunnel  diameter  for  two  different  strength  clays.  The 
Chicago  data  is  for  a tunnel  constructed  in  1940  where  the  soft  clay  had 
a shear  strength  of  700  psf  (34  kN/m2),  Terzaghi  (1943).  The  London 
data  is  for  a tunnel  constructed  in  1957  where  the  stiff  clay  had  a 
shear  strength  of  about  8 Ksf  (380  kN/m2),  Ward  and  Thomas  (1965). 
Although  there  is  much  scatter  in  this  type  of  data,  the  general  trend 
is  usually  the  same.  A steady  increase  occurs  with  sudden  changes  in 
the  observations  caused  by  a change  in  air  pressure  or  a passing 
ne i ghbor  tunne 1 . 

Of  course,  the  trend  of  liner  loads  increasing  with  time  must 
logically  halt  at  some  time  if  failures  are  to  be  prevented.  This  end 
point  when  liner  loads  become  constant  is  not  known  for  any  of  the 
reported  data,  since  the  measurements,  were  not  extended  over  a long 
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Figure  2.8;  Variation  of  Horizontal  Diameter  with  Time 
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enough  period  of  time.  If  this  effect  were  due  entirely  to 
consolidation,  it  should  be  expected  that  the  load  and  diameter  changes 
would  cease  when  the  consolidation  process  is  completed.  However,  there 
are  likely  several  causes,  such  as  creep,  leakage  into  the  liner  in 
addition  to  consolidation  and  the  problem  is  not  clearly  understood. 

2.4.3  Summary 

Field  data  documenting  the  effect  of  pore  pressure  dissipation  on 
tunnel  behavior  is  scarce.  The  data  available  suggests  there  are 
significant  ground  movements  and  liner  load  changes  caused  by 
consolidation  effects.  However,  at  this  point  the  cause  of  this 
behavior  is  only  poorly  understood.  In  fact,  there  are  a number  of 
factors  that  likely  contribute  to  the  time-dependent  behavior  of  tunnels 
in  clay.  It  would  appear  that  these  effects  may  be  examined  using  the 
finite  element  method. 


2 . 5 FINITE  ELEMENT  MODELLING  OF  THE  CONSTRUCTION  PROCEHURE 

The  finite  element  method  is  a very  useful  tool  for  analyzing  stress 
conditions  in  a body  of  soil  which  may  have  varying  soil  properties  and 
be  subjected  to  complex  boundary  conditions.  E 1 asto-p 1 ast i c material 
laws  may  be  used,  pore  pressures  may  be  calculated  along  with  stresses, 
strains  and  displacements.  The  tunnel  excavation  and  liner  installation 
procedure  may  be  simulated.  With  calibrations  against  field  behaviors, 
finite  element  programs  using  e 1 asto-p 1 ast i c parameters  can  be  used  for 
reliable  prediction  of  ground  movement  problems. 
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2.5.1 


Plane  Strain  Assumption 


A tunnel  heading  is  most  accurately  modelled  using  a 
three-dimensional  mesh.  However,  this  becomes  extremely  costly  for  an 
e 1 asto-p 1 ast i c consolidation  analysis.  It  is  possible  to  reduce  the 
problem  to  a two-dimensional  problem  under  certain  conditions. 

First,  it  is  assumed  that  movements  into  the  face  of  the  tunnel  are 
controlled.  This  can  readily  be  done  using  an  advanced  shield  or 
compressed  air.  Second,  it  is  assumed  that  the  permeability  of  the  soil 
is  such  that  the  construction  procedure  takes  place  in  undrained 
conditions.  Thus,  the  tunnel  heading  will  have  advanced  a considerable 
distance  before  the  excess  pore  pressures  begin  to  dissipate,  so  that 
the  consolidation  process  will  take  place  under  plane  strain  conditions. 
This  assumption  is  quite  reasonable  since  advances  of  65  feet  (20  m)  per 
day  can  be  achieved  and  for  most  clays  a negligible  amount  of 
consolidation  will  occur  in  the  first  few  days. 

2.5.2  Tunnel  Excavation 

First,  initial  stresses  are  set  up  in  the  finite  element  mesh,  (see 
Figure  2.9).  Then  the  excavation  of  the  tunnel  is  simulated  by 
relieving  the  in-situ  stresses  around  the  tunnel  opening.  Equivalent 
loads  are  applied  at  the  nodal  points  surrounding  the  opening  in  order 
to  model  this  excavation.  A detailed  account  of  this  procedure  is  given 
by  Mana  (1978). 
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2.5.3  L i ner  I nstal 1 at i on 


The  liner  is  modelled  using  eight  noded  isoparametric  elements  which 
have  been  found  to  give  good  results  when  subjected  to  bending  forces. 
The  stiffness  of  the  liner  is  reduced  to  account  for  the  flexibility 
introduced  by  using  bolted  liner  segments  rather  than  a cast  in  place 
concrete  liner.  The  liner  elements  are  in  place  from  the  start  of  the 
analysis,  but  are  initially  inactive.  They  are  activated  when  the  nodal 
point  at  the  crown  of  the  tunnel  has  moved  sufficiently  to  simulate  a 
closure  of  the  tail  void.  Using  this  procedure  the  size  of  the  tail 
void  may  be  altered  to  simulate  shield  design  and  construction 
performance.  The  weight  of  the  liner  is  included  in  the  analysis,  since 
the  weight  of  a concrete  liner  is  equal  to  a significant  proportion  of 
the  weight  of  the  excavated  soil. 

2.5.4  Soil  Model 

An  e 1 asto-p 1 ast i c soil  model  is  needed  to  provide  a good 
representation  of  the  stress-strain  behavior  of  clay,  since  it  will 
undergo  yielding  as  it  moves  into  the  tail  void.  The  Cam  clay  model, 
which  has  been  chosen  for  this  work,  has  been  developed  at  Cambridge 
University  over  the  past  15  years.  This  model  uses  a capped  yield 
surface  and  a hardening  rule  based  on  isotropic  compression  theory. 
Gunn  (1977)  has  used  Cam  clay  in  a finite  element  code  to  predict  the 
drained  behavior  of  tunnels  in  clay.  His  results  were  compared  with 
model  test  results  performed  by  Orr  (1976)  and  Seneviratne  (1977). 
Agreement  was  observed  between  the  finite  element  and  test  results  for 
surface  settlements  and  for  contours  of  volumetric  and  shear  strain. 
Details  of  the  Cam  clay  model  used  are  given  in  Chapter  III. 
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2.5.5  Conso 1 i dat i o n 


Excess  pore  pressures  are  built  up  during  excavation  and  liner 
installation.  The  construction  is  assumed  to  take  place  in  undrained 
conditions.  The  excess  pore  pressures  are  then  allowed  to  dissipate. 
The  surface  settlements  and  liner  stresses  are  altered  as  a result  of 
this  consolidation.  A detailed  discussion  of  the  consolidation  theory 
used  is  presented  in  Chapter  IV. 

2.5.6  Summary 

The  construction  of  a tunnel  in  clay  may  be  modelled  in  plane 
strain.  The  in-situ  stresses  are  setup  first,  and  this  is  followed  by 
excavation  of  the  soil.  The  tunnel  liner  is  then  activated  when 
movement  of  the  soil  into  the  tunnel  opening  has  closed  the  tail  void. 
The  finite  element  program  uses  an  e 1 asto-p 1 ast i c soil  model.  Excess 
pore  pressures  are  calculated  by  adding  pore  pressure  degrees  of  freedom 
to  the  nodal  points.  Loads  are  applied  in  small  increments,  and  the 
solution  remains  accurate  when  strains  become  large. 
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Figure  2.9:  Procedure  for  Simulation  of  Tunnel  Construction 
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Chapter  3 


DEVELOPMENT  OE  THE  CAM  CLAY  SOIL  MODEL  EOR  USE  IN  A FINITE  ELEMENT  CODE 
3 . 1 INTRODUCTION 

In  this  chapter  a path  will  be  traced  from  the  equations  that  define 
the  Cam  clay  soil  model  to  the  development  of  FORTRAN  subroutines  that 
calculate  the  stress-strain  matrix  for  use  in  a finite  element  code. 

The  finite  element  code  is  designed  to  obtain  solutions  to  problems 
of  combined  non-linear  material  and  geometric  behavior,  in  particular, 
problems  related  to  the  construction  of  tunnels  in  saturated  clays.  To 
solve  the  non-linear  equilibrium  equations,  incremental  element 
stiffness  matrices,  [K],  are  calculated  and  assembled  into  a global 
incremental  stiffness  matrix.  The  corresponding  equations  are  solved  to 
determine  the  incremental  di spacements.  Element  stiffness  matrices  are 
calculated  from 


[ K 1 


[ B J T [ D ] [ B ] d V 


V 


(3.  1) 


where  [B]  is  a matrix  relating  strains  to  displacements, 

[ D ] is  a matrix  relating  stresses  to  strains. 

Incremental  stresses  may  be  calculated  from  these  displacements  and 
added  to  the  stresses  that  existed  at  the  start  of  the  increment,  Marcal 
(1969),  Zienkiewicz  (1977). 
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Thus,  to  calculate  the  incremental  stiffness  matrix,  it  is  necessary 
to  determine  the  relationship  between  incremental  stresses  and 
incremental  strains 

{da'}  = [DHdd  (3.2) 

where  {da'}  is  a vector  of  incremental  effective  stresses, 

{de}  is  a vector  of  incremental  strains. 

In  accordance  with  the  standard  sign  convention  for  soil  mechanics, 
compressive  stresses  are  taken  as  positive  and  a reduction  in  volume  is 
taken  as  a positive  volumetric  strain. 

In  a cartesian  coordinate  system  the  effective  stress  vector  is 
given  by 

{a'}  = (ax'  ay'  az'  ryz  tzx  rxy)T  (3.3) 

The  corresponding  strain  vector  is  given  by 

{c}  = (ex  € y e z yyz  7zX  yxy)T  (3.4) 

The  pore  water  pressure  does  not  affect  the  effective  stress-strain 
matrix.  The  total  stress  is  given  in  the  usual  way  for  direct  stresses 
as 


a = a'  + u 


where  u is  the  pore  water  pressure. 


(3.  5) 
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3 • 2 THE  CAM  CLAY  MO  DEL 


3.2.1  I ntroduct  i on 

To  predict  the  response  of  soil  it  is  necessary  to  use  a 
mathematical  model.  The  earliest  researchers  chose  an  elastic  model. 
Soil  was  first  modelled  as  an  e 1 asto-p  1 ast i c work  hardening  model  by 
Drucker,  Gibson  and  Henkel  (1957).  Using  this  idea  in  conjunction  with 
critical  state  theory,  the  Cam  clay  model  was  developed  by  Roscoe  and 
Schofield  (1963),  and  by  Roscoe,  Schofield  and  Thurairajah  (1963).  The 
yield  surface  was  later  changed  from  the  original  bullet  shape  to  an 
elliptical  shape  in  Modified  Cam  clay,  Roscoe  and  Burland  (1968).  This 
is  the  version  of  the  Cam  clay  model  that  will  be  used. 

The  state  of  a point  in  a soil  continuum  is  defined  by  two 
generalized  stress  parameters,  p'  and  q,  and  the  voids  ratio  e,  Gunn 
(1977).  The  mean  effective  stress  is  defined  as 

<yx'  + ay'  + crz' 

p'  = (3.6) 

3 

The  octahedral  shear  stress  is  defined  by 
1 

q = — \Z(ay/-a2/)2+(crz,-CTx/)2+(ffx/-<jy,)2  + 6TyZ2  + 6Tzx2  + 6TXy2  (3.7) 

1/2 


3.2.2  Yield  Surface 

If  an  element  of  soil  is  loaded  until  it  reaches  the  yield  surface, 
any  further  loading  will  cause  yielding  so  that  the  modulus  of  the  soil 
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will  change  as  the  yield  surface  expands. 


If  this  element  is  then 


unloaded  the  soil  will  behave  elastically  again.  This  yielding  does  not 
coincide  with  failure,  which  occurs  when  the  soil  is  unable  to  take  any 
further  load. 

The  e 1 asto-p 1 ast i c material  behaves  elastically  provided  its  state 
in  p'-q-e  space  remains  below  the  yield  surface.  The  projection  of  this 
surface  on  the  p'-q  plane  for  a given  preconsolidation  pressure  p'c  is 
an  ellipse.  Many  other  shapes  of  the  yield  surface  in  the  p'-q  plane 
have  been  proposed.  The  earliest  work  by  Drucker,  Gibson  and  Henkel 
(1957)  used  a spherical  cap.  Early  work  in  Cam  clay  used  a bullet 
shape,  Schofield  and  Wroth  (1968).  An  exponential  cap  was  used  by 
DiMaggio  et  al  (1971),  and  a conical  cap  has  been  developed  more 
recently  by  Prevost  and  Hoeg  (1975).  Figure  3.1  shows  the  yield  curve 
whose  equation  is 

92 

F = — + p'(p'-p'c)  = 0 (3.8) 

M2 

where  M is  the  slope  of  the  failure  line.  The  significance  of  this  line 
is  explained  in  the  next  section. 

3.2.3  Fai lure  Enve 1 ope 

Failure  may  be  said  to  occur  when  strains  continue  to  increase 
without  any  increase  in  stresses.  This  happens  when  yielding  coincides 
with  the  critical  state  so  that 

q = Mp ' (3.9) 
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Figure  3.1  also  shows  the  failure  line. 


When  q < Mp'  stable  strain 


hardening  occurs.  Thus,  if  an  outward  probe  is  made,  the  yield  surface 
will  expand.  When  q > Mp ' unstable  strain  softening  occurs.  In  this 
case,  if  an  outward  probe  is  made,  the  yield  surface  collapses. 

The  finite  element  code  was  developed  for  making  predictions  in 
problems  for  which  q < Mp ' and  thus  produces  positive  definite  stiffness 
matrices.  Critical  state  may  also  be  handled  by  the  program  provided 
that  the  problem  is  strain  controlled.  Thus,  stress  controlled  problems 
that  reach  critical  state  and  strain  softening  problems  were  not  solved 
in  this  work. 

3.2.4  Flow  Rule 

The  normality  principle  is  used  so  that  the  plastic  strain  increment 
vector  is  normal  to  the  yield  surface.  This  is  always  the  case  for 
associative  plasticity.  Hill  (1950).  Recently  it  has  been  suggested, 
Banerjee  and  Stipho  ( 1978),  that  a non-assoc  i at  i ve  law  should  be  used. 
This  recommendation  is  based  on  experimental  evidence  by  Lewin  and 
Burland  (1970),  who  found  that  the  plastic  strain  increment  vector 
deviates  from  the  outward  normal  to  the  yield  surface.  Nevertheless, 
for  many  cases  associative  plasticity  has  given  good  results  , and 
because  it  leads  to  a symmetric  stiffness  matrix,  the  cost  savings  are 
significant. 
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Shear  Stress, 


Figure  3.1:  Cam  clay  Yield  Model 
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3.2.5  Hardening  Rule 


When  yielding  occurs,  the  size  of  the  elliptical  yield  curve  is 
altered  since  pc',  the  hardening  parameter,  changes.  This  parameter  is 
governed  by  the  virgin  isotropic  consolidation  curve  shown  in  Figure  3.2 
and  is  given  by 

e=e1-AlogePc'  (3.10) 

where  e i is  the  value  of  the  voids  ratio  for  unit  p c ' , 

A is  the  slope  of  the  virgin  isotropic  consolidation  curve. 

In  practice  this  curve  may  deviate  slightly  from  a straight  line 
when  plotted  on  semi-log  graph  paper.  For  this  reason,  Chang  and  Duncan 
(1977)  chose  to  interpolate  the  line  from  a set  of  experimental  data  fed 
into  their  program. 

3 . 3 FLASTIC  EFFECTIVE  STRESS-STRAIN  BFHAV I OR 

Below  the  yield  surface,  Cam  clay  is  a non-linear  elastic  material. 
The  bulk  modulus,  K,  may  be  determined  from  the  elastic  rebound  line, 
(see  Figure  3.2),  and  is  given  by 

1+e 

K = p'  (3.11) 

K 

where  K is  the  slope  of  the  rebound  isotropic  consolidation  curve.  One 
further  elastic  property  needs  to  be  specified,  since  the  others  are 
related  by 

E = 3K( 1-2v)  = 2G( 1+v)  (3.12) 
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where  E is  Young's  modulus, 

G is  the  shear  modulus, 
v is  Poisson's  ratio. 


For  some  problems  it  may  be  convenient  just  to  specify  Poisson's 
ratio;  a value  of  0.3  for  clays  usually  gives  satisfactory  results. 
However,  it  is  preferred  that  the  shear  modulus  be  specified.  This 
value  need  not  be  a constant,  but  may  vary  with  depth.  This  method  was 
used  for  the  tunnel  analyses  in  this  report. 


Thus,  the  elastic  stress-strain  matrix  may  be  given  by,  Timoshenko 
and  Goodier  (1970), 


1 

[Dei  = - 
3 


3K+4G  3K-2G 
3K-2G  3K+4G 
3K-2G  3K-2G 
0 0 

0 0 

0 0 


3K-2G  0 

3K-2G  0 

3K+4G  0 

0 3G 

0 0 

0 0 


0 O' 

0 0 

0 0 

0 0 

3G  0 

0 3G- 


(3.13) 
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Mean  Effective  Stress,  p' 

Figure  3.2:  Cam  clay  Hardening  Model 
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3 . 4 ELASTO-PLASTIC  EFFECTIVE  STRESS-STRAIN  BEHAVIOR 


3.4.1  Derivation  of  the  Stress-Strai n Matrix 

It  is  now  necessary  to  derive  an  e 1 asto-p 1 ast i c stress-strain  matrix 
for  a capped  yield  model.  The  yield  surface  given  in  equation  (3.8)  may 
be  rewritten  as 

F(  {a'}T,Pc')  =0  (3.14) 

Taking  the  total  derivative  gives 


aT  <>F 

(do'}  + dpc'  = 0 

Jpc' 


The  flow  rule  is  stated  as 


{deP} 


(3.15) 


(3. 16) 


where  4>  is  a proportionality  constant,  and  superscripts  e and  p refer  to 
the  elastic  and  plastic  components  respectively.  Since  the  volumetric 
plastic  strain  increment  is  given  by 

devp  = dexP  + deyP  + dc2P  (3.17) 

The  flow  rule  may  be  rewritten  as 

bF 

de VP  - <f>  (3.18) 

*P' 

The  total  strain  is  the  sum  of  the  elastic  and  plastic  components, 
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{de}  = {de*}  + {deP } 


(3.19) 


Substituting  equations  (3.2)  and  (3.16)  into  (3.19)  gives 


{de } = [ De ] ' 1 {da' } + * 


-bcr 


(3. 20) 


From  (3. 15),  (3. 18)  and  (3.20) 


bF  T T 


l De  ] 


-bo 


{de } - 4> 


ba'JJ 


bF  dp c'  bF 

+ <t>  — = 0 

bpc'  devP  2>p' 


(3.21) 


Solving  equation  (3.21)  for  $ and  substituting  into  equation  (3.20) 
g i ves 


Iba'J  Iba'J 


l D ] = { ] - 


(3.22) 


rbF  TT  rbF 

[ De  ] 

-ba'J  tba 


bF  dp c'  bF 


bp c ' devP  bp' 


This  is  the  final  form  of  the  e 1 asto-p 1 ast i c stress-strain  matrix^ 


3.4.2  Calculation  of  the  Derivatives 

Equation  (3.22)  is  a general  form  of  the  incremental  effective 
stress-strain  matrix  for  any  e 1 asto-pl ast i c material  with  a capped  yield 
model.  To  evaluate  this  matrix  for  the  Cam  clay  model  it  is  necessary 
to  calculate  all  of  the  derivatives  in  this  equation. 

To  evaluate  the  vector  {bF/bo'},  the  chain  rule  for  partial 
differentiation  may  be  used  as  follows 
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Differentiating  equation  (3.5)  with  respect  to  o'  gives 


' 1/3' 
1/3 
1/3 
0 
0 

. 0 . 


Differentiating  equation  (3.7)  with  respect  to  a'  gives 


f 3(a/-pr)' 

[?L1 

3(<xy'-p') 
3(ff z/-P' ) 

Us'. 

2q 

6 T y 2 

6tzx 

■ 6^  x y 

Differentiating  the  yield  curve,  equation  (3.8),  with  respect 
gives 

bF 

= 2p'  - Pc' 

2>P' 


(3.24) 


(3.25) 


to  p ' 


(3.26) 


Differentiating  the  yield  curve,  equation  (3.8),  with  respect  to  q gives 


Substituting  the  terms  calculated  in  (3.24)  through  (3.27)  into  equation 


(3.23)  gi ves 


' 2p'-pc' 
3 

2p'-pc' 

3 

2p'-pc' 

3 


3(<rx'-p/)  ' 
M2 

3(ay,-p/) 

m2 

3(ciz,-p/) 

M2 


6T  yZ 


M2 


STj-x 

M2 

6t  Ky 
M2 


(3.  28) 


The  remaining  derivatives  to  be  calculated  concern  the  work 
hardening  parameter.  Differentiating  the  yield  curve,  equation  (3.8), 
with  respect  to  pc'  gives 


bF 

= -p' 

2>Pc' 


The  final  derivative  may  be  calculated  as  follows 


(3.29) 


dpc'  dpc'  de 

= (3.30) 

devp  de  devp 

Differentiating  equation  (3.11)  gives 
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dpc'  Pc' 

= (3.31) 

de  A 

And  the  relationship  between  voids  ratio  and  plastic  volumetric  strain 
i s 

de  A ( 1+e0 ) 

= (3.32) 

devp  A-k 


where  e<>  is  the  voids  ratio  at  the  start  of  the  increment.  Thus,  the 
hardening  parameter.  A,  is  given  by 


bF  dpc'  'dF 

A = (3.33) 

dpc'  devp  dp' 


or 


p/pc'(2p'-pc,)(1+e0) 

A = (3.34) 

A-* 

From  these  equations  the  stress-strain  matrix  may  be  calculated  for 
the  Cam  clay  model  when  the  yield  surface  has  been  reached. 

3. 5 MATERIAL  PARAMETERS  REQUIRED 

For  a soil  model  to  be  useful,  the  material  parameters  that  describe 
the  model  should  be  few  and  easily  obtainable  from  conventional 
laboratory  tests.  Cam  clay  may  be  described  by  the  following  five 
parameters : 
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1.  X,  the  slope  of  the  virgin  isotropic  consolidation  curve 


2.  K » the  slope  of  the  rebound  isotropic  consolidation  curve 


3.  ecs,  the  voids  ratio  for  unit  p'  on  the  critical  state  line 


4.  M,  the  slope  of  the  failure  line  in  p'-q  space 


5.  G,  the  elastic  shear  modulus 


The  first  three  parameters  may  be  readily  obtained  from  an  isotropic 

consolidation  test.  The  slope  of  the  failure  line,  M,  may  be  calculated 

from  the  friction  angle  , obtained  from  a set  of  drained  tri  axial 
tests,  since 

6 sin  4>' 

M = (3.35) 

3 - sin 

The  final  elastic  constant,  G,  for  a particular  confining  stress  03', 
can  be  determined  as  one  third  of  the  initial  gradient  of  the  deviatoric 

stress  axial  strain  curve  of  an  undrained  triaxial  test,  Wroth  and 

Zytynski  (1977). 

Since  the  shear  modulus  usually  varies  with  depth  in  practice,  using 
a single  value  for  this  parameter  is  a simplification.  The  shear 
modulus  may  be  written 

G = mcu  ( 3.  36) 

where  cu  is  the  shear  strength, 
m is  a multiplier. 
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The  multiplier  may  be  assumed  to  be  a constant.  The  shear  strength 
may  be  assumed  to  vary  linearly  with  depth;  for  example,  it  is  equal  to 
approximatel  y 0.3  times  the  effective  vertical  stress  for  many  normally 
consolidated  clays.  A further  sophistication  may  be  added  by  assuming 
that  a desiccated  layer  exists  above  any  normally  consolidated  zone, 
thus  prohibiting  the  value  of  G from  becoming  zero  at  the  surface.  This 
is  the  type  of  approach  used  in  the  tunnel  analyses. 

There  are  other  material  parameters  which  are  needed  in  order  to  set 
up  the  initial  stresses  in  a problem.  The  following  parameters  are 
required  for  this  purpose: 

1.  7,  the  density  of  the  soil 

2.  ko,  the  coefficient  of  earth  pressure  at  rest 

3.  OCR,  the  overconsolidation  ratio 


Details  as  to  how  the  initial  stresses  are  set  up  in  the  program  will 
be  given  in  the  chapter  VI. 

3.  6 REDUCTION  TO  ONE-  AND  TWO-DIMENSIONAL  PROBLEMS 

The  effective  stress-strain  matrix  has  been  derived  to  this  point  in 
the  chapter  in  its  most  general  form  for  use  in  three-dimensional 
problems.  However,  for  many  applications  the  problem  reduces  to  two 
dimensions;  for  example,  a long  embankment  may  be  analyzed  as  a plane 
strain  problem,  or  a triaxial  test  may  be  treated  as  an  axisymmetric 
problem.  Some  problems  even  reduce  to  a single  dimension;  for  example. 
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consolidation  of  a thin  clay  layer  under  a large  area  of  uniform  load. 


Other  problems  approximate  to  these  conditions,  such  as  consolidation 
around  a tunnel.  Thus,  for  the  sake  of  simplicity  and  considerably 
reduced  computer  costs,  it  is  often  desirable  to  solve  problems  of  fewer 
than  three  dimensions. 

3.6.1  Ax i symme  tr i c 

Consider  a problem  which  has  the  y axis  as  the  axis  of  symmetry. 
Then  ryz,  and  rzx  may  be  set  to  zero.  For  this  reason  the  original  6x6 

stress-strain  matrix  reduces  to  a 4x4. 

/ 

3.6.2  Plane  Strain 

For  a problem  with  the  longitudinal  axis  aligned  along  the  z axis, 
ryz  and  t z x may  be  set  to  2ero.  And  although  ez  is  zero,  the  value  of 
az  is  not  necessarily  zero.  Also  the  relationship 

<jz  = v(<rx+cxy)  ( 3.  37) 

does  not  hold  for  an  e 1 asto-p 1 ast i c material  such  as  Cam  clay.  Thus, 
this  problem  also  reduces  to  a 4x4. 

The  importance  of  cr2  for  analysis  of  deep  circular  tunnels  was  shown 
by  Pender  (1980).  He  found  that  pore  pressure  changes  during  the 
construction  of  a lined  tunnel  depended  on  cr  2 > and  differred  drastically 
from  those  predicted  using  the  Stress  Path  Method,  (Lambe  and  Marr 
1979),  where  a z is  neglected.  Early  versions  of  the  Cam  clay  model 
which  were  used  in  finite  element  programs  also  neglected  the  effect  of 
c 2 (Simpson  and  Wroth  1972).  The  effect  of  is  included  in  the  tunnel 
analyses  in  this  report. 
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3.6.3  One  Dimension 


All  the  shear  stresses  are  zero  in  a problem  that  reduces  to  one 
dimension.  Thus,  the  stress-strain  matrix  reduces  to  a 3x3.  Further 
simplifications  could  be  made  when  initial  stress  conditions  allow  two 
of  the  direct  stresses  to  be  identical. 

3.7  FULLY  DRAINEP  AND  UNDRAINEP  RFHAVIOR 

In  chapter  IV  the  soil  model  is  expanded  to  include  a consolidation 
capability  by  adding  a pore  pressure  degree  of  freedom  to  the  nodal 
points.  However,  fully  drained  and  undrained  problems  may  be  solved  in 
a more  simple  way.  The  equations  presented  in  this  chapter  apply  for  a 
drained  analysis.  For  an  undrained  analysis,  when  calculating  the 
stress-strain  matrix,  the  bulk  modulus  of  water  is  added  to  the  bulk 
modulus  of  the  soil.  This  bulk  modulus  is  a very  large  number  when 
compared  with  the  shear  modulus  of  the  soil.  Thus,  the  effect  is  to 
produce  no  volume  change  (for  numerical  reasons  a very  small  amount  of 
volume  change  is  necessary),  but  to  allow  shearing  to  take  place.  This 
method  is  preferable  to  the  usual  method  of  setting  Poisson's  ratio 
c 1 ose  to  0.5. 

During  the  stress  recovery  cycle,  the  stress-strain  matrix  should  be 
computed  without  adding  the  bulk  modulus  of  water.  The  effective  stress 
increment  is  then  obtained.  The  pore  pressure  may  be  obtained  by 
multiplying  the  small  volume  change  by  the  bulk  modulus  of  water. 
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3.8  FINITE  ELEMENT  ALGORITHMS 


3.8.1  Introduction 

FORTRAN  subroutines  for  calculating  the  stress-strain  matrix  for  Cam 
clay  and  for  calculating  the  updated  stresses  in  the  stress  recovery 
cycle  are  given  in  Appendix  B. 

It  is  useful  to  understand  this  part  of  the  program  design.  During 
the  assembly  of  the  global  stiffness  matrix,  stiffness  matrices  are 
calculated  for  each  element  using  the  numerical  technique  of  gaussian 
integration.  The  stress-strain  matrix,  [D],  is  calculated  at  each 
integration  point.  Since  this  matrix  depends  on  the  stress  state  and 
history  of  that  integration  point,  these  stresses  must  be  stored  from 
one  increment  to  the  next. 

3.8.2  Stress-Strain  Matrix 

The  stress-strain  matrix  subroutine  is  called  for  each  integration 
point  and  computed  according  to  equation  (3.22).  First  the  elastic 
component  , [Del,  is  computed.  The  stress  state  is  examined,  and  if  it 
is  below  the  yield  surface,  the  behavior  is  elastic  and  thus  no  further 
calculations  are  required.  If  the  stress  state  lies  on  the  yield 
surface,  the  remainder  of  the  formula  in  equation  (3.22)  must  be 
cal cul ated . 

For  cyclic  problems,  the  stress  state  may  lie  on  the  yield  surface 
even  though  the  soil  is  unloading.  To  accommodate  this  possibility,  it 
is  necessary  to  be  allow  the  program  to  repeat  or  iterate  an  increment 
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until  the  guess  (unloading  or  loading)  made  in  the  assembly  cylcle  is 


matched  by  the  correct  response  in  the  recovery  cycle.  Since  the 
problems  analyzed  in  this  work  are  generally  not  cyclic,  this  capability 
has  been  omitted.  The  problem  may  be  overcome,  where  necessary,  by 
taking  very  small  load  increments. 

3.8.3  Stress  Updating 

The  stress-strain  matrix  and  the  strain-displacement  matrix  are 
recovered  from  temporary  storage,  where  they  are  preserved  from  the 
assembly  cycle.  The  incremental  stresses  are  then  calculated  from 

{do'}  = iDllBHdx}  (3.38) 

where  { d x } is  a vector  of  the  incremental  displacements.  These 
incremental  stresses  are  added  to  the  total  stresses  which  existed  at 
the  start  of  the  increment.  The  yield  surface  is  updated  if  it  has 
expanded  during  the  increment.  The  updated  voids  ratio  is  calculated 
from 


e = e cs  - K 1 og e p' 


(X-x) 


(3.39) 


The  updated  quantities  p',  q,  e and  pc'  are  saved  for  each  integration 
point  along  with  the  stresses  and  strains. 


3.9  COMPARISON  OF  THE  FINITE  ELEMENT  RESULTS  WITH  TEST  DATA 

At  this  stage  it  is  desirable  to  test  the  response  of  the  finite 
element  code  using  the  Cam  clay  model.  This  is  done  by  analyzing  the 
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problem  of  drained  triaxial  tests  on  normally  consolidated  Weald  clay. 


The  experimental  tests  were  performed  by  Bishop  and  Henkel  (1971).  The 
material  properties  used  are  shown  in  Table  3.1. 


Weald  Clay 

TABLE  3. 1 

Material  Parameters 

K 

0.031 

A 

0.088 

£ C S 

1 .0575 

M 

0.882 

G 

434.8  psi 

The  finite  element  mesh  used  is  shown  in  Figure  3.3.  The  test  was 
simulated  using  one  four  noded  axisymmetric  isoparametric  element.  The 
initial  stresses  were  set  equal  to  the  value  of  03',  which  was  30  psi, 
and  then  the  test  was  run  by  applying  displacement  increments  in  the 
negative  y direction  to  nodes  1 and  2.  Thus,  a strain  controlled  test 
was  simulated,  and  hence  problems  of  instabilty  at  failure  were  avoided. 
Results  are  shown  in  Figure  3.4  and  excellent  agreement  was  obtained 
between  experimental  and  finite  element  results. 
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Figure  3.3:  Finite  Element  Mesh  for  Triaxial  Tests 
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Figure  3.4:  Drained  Triaxial  Compression  Tests  on  Weald  Clay 
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Chapter  4 


EXTENDING  THE  MODEL  TO  INCLUDE  CONSOLIDATION 
4.  1 INTRODUCTION 

The  settlement  of  soils  under  load  due  to  the  squeezing  out  of  the 
pore  water  is  a process  called  consolidation.  This  phenomenon  was  first 
explained  by  Terzaghi  (1925)  for  the  case  of  one-dimensional 
consolidation.  The  theory  was  later  generalized  to  three  dimensions  by 
Biot  (1941). 

A solution  to  the  problem  of  consolidation  involves  coupling  the 
equations  of  elasticity  with  those  of  fluid  flow.  To  solve  these 
equations  using  a finite  element  scheme,  the  problem  must  be  discretized 
in  both  the  space  and  time  domains.  Sandhu  and  Wilson  (1969)  proposed  a 
finite  element  procedure  for  the  case  of  seepage  in  an  elastic  medium. 
Small,  Booker  and  Davis  (1976)  used  a different  technique  to  find  a 
solution  to  the  problem  of  e 1 as to-p 1 as t i c consolidation  of  soil.  It  is 
this  method  that  is  used  here,  since  Cam  clay  is  a work  hardening 
e 1 asto-p 1 ast i c material.  Only  the  essential  parts  of  the  theory  will  be 
outlined;  readers  interested  in  the  complete  mathematical  solution  are 
referred  to  appropriate  papers  in  the  text. 
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4 . 2 GOVERNING  EQUATIONS 


In  this  section  the  partial  differential  equations  governing  the 
e 1 asto-p  1 ast  i c consolidation  of  soil  will  be  presented.  Appropriate 
boundary  conditions  will  also  be  given.  Double  suffix  notation  is  used, 
and  a superior  dot  denotes  material  differentiation  with  respect  to 
time. 

The  governing  equations  are 

1.  Equilibrium  between  the  incremental  stresses  and  the  incremental 
body  forces: 

Sa  i j 

F ? = 0 (4.1) 

dx  j 

where  a , j are  the  components  of  the  total  stress  tensor, 

Xj  are  the  coordinates  in  a cartesian  reference  frame, 

Fi  are  the  components  of  body  force. 

2.  Stress-strain  relation  between  the  incremental  effective  stresses 
and  the  incremental  strains: 

" pSij  = Djjkl€kl  (4.2) 

where  cn  j ' are  the  components  of  the  effective  stress  tensor, 
p is  the  pore  water  pressure, 

S,j  is  the  Kronecker  delta, 

D,-jki  are  the  coefficients  of  the  stress-strain  law,, 
eki  are  the  components  of  the  strain  tensor. 
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3.  Darcy's  law  for  the  rate  of  flow  of  water  through  soil: 


kij  Sp 

v*  (4.3) 

y w Sx  3 

where  V3  are  the  components  of  the  superficial  velocity  vector, 
kij  are  the  coefficients  of  the  permeability  matrix, 
y w is  the  unit  weight  of  water. 

4.  And  finally,  relationship  between  the  rate  of  volume  decrease  of 
the  element  and  the  rate  at  which  water  is  expelled,  assuming  the 
pore  water  is  incompressible  relative  to  the  soil  skeleton: 

bv  3 

= (4.4) 

2>x  3 M 

where  €v  is  the  volumetric  strain, 
t is  the  time. 


To  model  problems  accurately,  many  different  types  of  boundary 
conditions  must  be  applied.  The  following  four  are  common  in  practice, 
and  may  be  used  in  the  finite  element  program: 

1.  To  simulate  a constrained  boundary,  the  displacement  increment  is 
set  to  zero: 

U3  = 0 (4.5) 

where  U3  are  the  displacements. 
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2.  To  allow  free  displacements,  the  applied  tractions  must  be  in 
equilibrium  with  the  stresses: 

(4.6) 

where  nj  is  the  outward  normal  to  the  surface, 

T;  are  the  applied  tractions. 

3.  To  simulate  a pervious  boundary,  the  pore  pressure  must  be 
allowed  to  remain  constant: 

P = 0 (4.7) 

4.  To  simulate  an  impervious  boundary,  water  must  not  be  allowed  to 
cross  the  boundary: 

n fv  * = 0 (4.8) 


A solution  to  these  equations  was  found  by  Booker  (1973). 


4 . 3 FINITE  ELEMENT  FORMULATION 

The  consolidation  problem  may  be  discretized  using  finite  elements. 
A pore  pressure  degree  of  freedom  is  added  to  the  usual  displacement 
degrees  of  freedom  at  the  nodal  points.  The  spatial  integrations  may 
then  be  performed  using  gaussian  quadrature,  and  a marching  process  is 
used  to  perform  the  time  integration. 
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Shape  functions  are  used  to  relate  the  continuous  value  of  the 
displacements,  u,  to  the  nodal  values,  { S } . 

(u)  = {NU}T {6}  (4.9) 

where  {Nu}  are  the  shape  functions. 

Shape  functions  may  also  be  used  to  relate  the  continuous  values  of 
the  pore  pressure,  p,  to  the  nodal  values,  (q). 

P = {Np}T{q}  (4.10) 

where  {Np}  are  the  shape  functions. 

Both  sets  of  shape  functions  depend  on  the  particular  type  of  finite 
element  used.  The  shape  functions  may  use  the  same  or  different  orders 
of  interpolation,  and  thus  a pore  pressure  degree  of  freedom  need  not 
necessarily  exist  at  each  nodal  point.  This  will  be  discussed  further 
i n Sect i on  4.10. 

The  coupled  equations  for  an  element  may  be  written  in  matrix  form, 
Carter,  Booker  and  Small  (1979): 


' K 

-LT  ' 

' dS  ' 

f 

- -L 

-Bdt$  - 

• dq  - 

« ndt  + $q0dt  - 

(4.11) 


where  the  element  stiffness  matrix,  [K],  is  given  in  the  usual  way  by 


[ K 1 


[ B IT [ D 1 [ B IdV 


V 


(4. 12) 
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The  matrix  [ L J couples  the  stiffness  equations  with  those  of  fluid  flow, 
and  is  given  by 


I L ] T 


{ B v) T (NP}TdV 


V 


(4.13) 


where  the  volumetric  strain-displacement  vector,  { B v } , is  given  by 


€ v = {Bv} {S } (4. 14) 

This  vector  is  easily  calculated  from  the  strain-displacement  matrix, 
since 

{Bv}  = ( 1 1 1 0 0 0 ) [B]  (4.15) 


The  fluid  flow  matrix  is  given  by 


1*1 


1 

y w 


[E]T[k]T[EldV 

V 

V 


where  [k]  is  the  permeability  matrix,  and  where 


I E IT  = 


(4. 16) 


(4. 17) 


Also,  { f } is  the  incremental  load  vector, 

{ q o } are  the  pore  pressures  at  the  beginning  of  the  increment, 

B is  an  integration  parameter,  which  is  discussed  in  more  detail 
i n Sect i on  4.11. 
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Finally,  the  existence  of  a hydrostatic  pore  pressure  distribution  is 


taken  care  of  by  adding  the  following  term  to  the  load  vector: 


Cn} 


[E]T[k]T{ig}dV 


V 


(4.18) 


where  {ig}  is  a vector  containing  the  direction  of  the  gravitational 

field. 


4.4  MATERIAL  PARAMETERS  REQUIRED 

In  order  to  extend  the  Cam  clay  model  to  include  the  effect  of 
consolidation,  only  two  further  material  properties  need  to  be 
specified.  These  are  as  follows: 

1.  y w , the  density  of  water  or  other  pore  fluid. 

2.  [k],  the  components  of  the  permeability  matrix. 

Since  the  density  of  water  is  known,  it  just  remains  to  find  the 
components  of  the  permeability  matrix.  These  may  be  obtained  from 
laboratory  or  field  tests.  Different  values  of  horizontal  and  vertical 
permeabilities  may  be  used. 


4 . 5 REDUCTION  TO  ONE-  AND  TWO-DIMENSIONAL  PROBLEMS 

In  Section  3.6  it  was  shown  that  for  problems  with  fewer  than  three 
dimensions,  the  computation  of  the  stress-strain  matrix  is  simplified, 
as  is  the  computation  of  the  fluid  flow  matrix,  [§].  Computation  of  the 
coupling  matrix,  [ L ] , however,  remains  unchanged. 
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For  problems  that  reduce  to  two  dimensions,  i.e.,  plane  strain  or 
ax i symmetr i c , the  order  of  [ E ] is  reduced,  since 


= 0 (4. 19) 

The  permeability  matrix  is  also  reduced,  so  that  for  horizontal  or 
vertical  bedding, 


Ik]  = 


(4.20) 


For  problems  that  reduce  to  one  dimension,  further  simplification  is 
permitted: 


[ E ]T 


(4.21) 


and 

[k ] = t k x i (4.22) 

4.6  FINITE  ELEMENT  ALGORITHMS 

A marching  process  is  used  to  perform  the  integration  through  time. 
This  involves  calculating  the  stress  state  at  some  future  time,  t+At, 
from  the  current  stress  state  at  time,  t,  due  to  the  excess  pore 
pressures  existing  at  time  t.  The  process  can  then  be  marched  forward 
by  taking  successive  time  increments.  At.  These  time  increments  are 
read  in  with  the  load  increment  data  for  each  increment. 
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At  each  increment  the  coupled  set  of  equations  needs  to  be 


recalculated.  This  is  done  using  gaussian  integration  over  each 
element.  The  element  matrices  are  then  assembled  into  a global 
stiffness  matrix.  A schematic  diagram  of  the  element  matrices  that  need 
to  be  calculated  is  shown  in  Figure  4.1.  It  should  be  noted  that  the 
element  stiffness  matrix  is  symmetric.  The  parts  A1,  A2  and  A3  are 
calculated  regardless  of  whether  the  consolidation  option  is  being  used. 
The  element  stiffness  matrix  is  stored  in  the  computer  program  in  a 
vector.  This  vector  stores  the  elements  of  the  stiffness  matrix  column 
by  column,  but  only  down  to  and  including  the  diagonal.  This  procedure 
has  two  advantages:  first  no  storage  is  wasted;  and  second,  when 
consolidation  is  included,  it  is  no  trouble  to  add  the  additional 
columns.  The  parts  Bl,  B2  and  B3  correspond  to  the  fluid  flow 
equations,  while  Cl  is  the  coupling  matrix.  FORTRAN  subroutines  for  the 
calculation  of  parts  Bl,  Cl  and  B3  are  given  in  Appendix  B. 


As  an  example,  the  matrices  needed  for  the  addition  of  pore 
pressures  will  be  shown  for  a simple  element.  The  element  chosen  is  an 
isoparametric  one-dimensional  element  with  three  nodes.  The  element  is 
named  L3P2  since  it  is  a link  element  with  displacement  degrees  of 
freedom  at  three  nodes,  and  pore  pressure  degrees  of  freedom  at  two 
nodes.  The  element  is  shown  in  Figure  4.2. 


The  shape  functions  that  describe  the  displacement  field  are 


{Nu> 


-£(1-£)/2  ' 
f ( 1 + £ )/2 


l 1-f* 

The  shape  functions  that  describe  the  pore  pressure  field  are 


(4.23) 


59 


A3 


= 1 — 

\ 1 

\ 

1 

\ 1 

\ ai  ! ci 

\ | 

\ 

\ 1 
\ 1 
\ | 

A2 

\ 

\ 1 

' \ 

| \ 

, \ 

\ 

\ B1 

| 

i \ 

i 

B2 

' 

1 \ 

1 <=J 

B3 


Figure  4.1:  Schematic  of  the  Element  Equations 
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Figure  4.2:  Element  L3P2 
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( 1-f )/2 


(4.24) 


{Np} 


(1+f)/2  J 


The  permeability  matrix  reduces  to  a 1x1  with  this  component  being  k*. 
If  there  is  no  acting  gravitational  field,  {ig}  is  zero.  The  flow 
matrix  becomes 


m 


+ 1 


k x 

■yw 


l E ]T  [ E ] |JP|  df 

V 

-1 


(4.25) 


where  |jp|  is  the  determinant  of  the  Jacobian  for  the  pore  pressure 
shape  functions,  and  is  given  by 


|JpI 


1 1 

xi  + -x2 

2 2 


(4.26) 


Also, 


[E]T 


-1/2  ' 
1/2  - 


The  coupling  matrix  becomes 


(4.27) 


+ 1 


( L ] T = 


[BV]T[NP]T  I Ju|  df 


-1 


(4.28) 


where  |ju|  is  the  determinant  of  the  Jacobian  for  the  displacement  shape 
functions,  and  is  given  by 
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*2  " 2 f X 3 


(4.29) 


|Ju|  = 

V 

s - - 

Xl  + 

1' 

f + ~ 

l 2-1 

l 2 i 

Since  there  are  two  Jacobians  for  this  one  element,  some  confusion  may 
arise  over  which  one  to  use.  When  node  3 is  placed  midway  between  nodes 
1 and  2,  the  values  of  both  Jacobians  become  identical. 


4 . 7 VERIFICATION  OF  ONE-DIMENSIONAL  SEEPAGE 

In  this  section  the  seepage  capability  of  the  program  is  tested. 
The  porous  medium  is  assumed  rigid  so  that  only  the  flow  character i st i cs 
are  tested.  Thus,  consolidation  does  not  take  place,  and  the  coupling 
matrix,  [ L 1 , is  not  tested.  Two  problems  are  analyzed;  the  first  is 
one-dimensional  seepage  and  the  second  is  radial  seepage. 


4,7.1  Vertical  Drainage 

A one-dimensional  flow  problem  was  analyzed.  The  pore  pressure  at 
depth  x = 10ft  is  given  as  p=1000psf,  and  at  depth  x = 1 0 0 f t it  is  given  as 
zero.  The  governing  differential  equation  is 


d 

dx 


k 


x 


dp' 

dx- 


= 0 


(4.30) 


Assuming  a constant  value  for  the  permeability,  k*,  the  equation  may 
be  integrated  twice  and  the  boundary  conditions  applied.  The  expected 
linear  distribution  of  pore  pressure  is  found  to  be 

100 

P = ( 100  - x ) (4.31) 
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The  finite  element  mesh  consisted  of  ten  L3P2  elements. 


Fourteen 


time  increments  were  used  to  allow  the  solution  to  stabilize.  The  size 
of  the  time  increment  increases  in  a logarithmic  fashion,  and  is 
dependent  on  the  value  of  kx  chosen.  The  results  are  shown  in  Figure 
4.3.  No  difference  is  visible  between  the  theoretical  results  arid  those 
from  the  finite  element  analysis. 


4.7.2  Radial  Drainage 

An  axisymmetric  flow  problem  was  analyzed.  The  pore  pressure  on  the 
inside  of  a t la  i c k porous  pipe  with  inside  radius,  r = 1 0 i n , is  p = 1 0 0 0 p s i , 
and  the  pressure  on  the  outside,  r=100in,  is  zero.  The  problem  is 
axisymmetric  and  also  plane  strain.  The  governing  differential 
equation,  written  in  polar  coordinates,  is 


d 

dr 


k 


r 


dp" 
r — 
dr- 


0 


(4.32) 


After  integrating  this  equation  twice  and  applying  the  boundary 
conditions,  the  pore  pressure  distribution  is  given  by 


P = 1000 


loger  ' 

2 

1 og e 1 0 - 


(4.33) 


The  finite  element  mesh  used  was  similar  to  that  used  in  the 
previous  problem,  The  same  time  increments  were  also  utilized.  The 
results  are  shown  in  Figure  4.4,  and  the  results  of  the  finite  element 
analysis  appear  to  agree  exactly  with  the  theory. 
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Pore  Pressure,  p (psf) 


Figure  4.3:  Pore  Pressure  Distribution  for  Vertical  Drainage 
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Pore  Pressure,  p (psf) 


Figure  4.4:  Pore  Pressure  Distribution  for  Radial  Drainage 
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4.8  V ERIFI  CATION  OF  ON E- D I MEN S I ONAL  CONSOLIDATION 


One-dimensional  consolidation  occurs  uhen  a uniform  load  is  applied 
to  a clay  layer  and  drainage  is  allowed  to  occur  in  the  vertical 
direction  only.  Thus,  settlement  occurs  with  time  as  the  excess  pore 
pressure  is  allowed  to  dissipate.  If  the  soil  skeleton  is  assumed  to 
behave  elastically,  the  consolidation  process  obeys  Terzaghi's 
consolidation  equation,  as  defined  by  Lambe  and  Whitman  (1969): 


(i^Ug  g b(7  z 

bz2  bt  bt 


(4.34) 


Thus,  the  excess  pore  pressure,  ue>  is  dependent  on  the  time,  t,  and  the 
depth,  2.  The  applied  vertical  stress  is  az,  and  cv  is  given  by 


k z 

cv  = — Ec  (4.35) 


where  Ec  is  the  constrained  Young's  modulus,  equal  to  E uhen  Poisson's 
ratio  is  zero.  The  solution  may  be  written  in  non-dimensional  form, 
Taylor  (1948),  so  that  the  excess  pore  pressure  is  related  to  the 
initial  pore  pressure,  uo»  by 

w 2u0  -M2T 

ue  = X sin(MZ)  e (4.36) 

m = 0 M 


where 


TT 

M = - ( 2m+ 1 ) (4.37) 
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The  non-dimensional  height  is  Z = z/H,  where  H is  the  drainage  distance. 
The  non-dimensional  time  is  given  by 

c vt 

T = (4.38) 

H2 

The  percentage  of  consolidation  for  the  clay  layer  is  given  by 
co  2 -M2T 

U = 1 - l — e (4.39) 

m = 0 M2 

The  finite  element  mesh  used  to  verify  the  consolidation  capability 
of  the  program  consisted  of  five  L3P2  elements,  and  is  shown  in  Figure 
4.5.  The  boundary  conditions  chosen  correspond  to  single  drainage, 
since  the  bottom  node  was  made  impervious  and  the  top  node  pervious. 
The  bottom  node  is  fixed  and  a uniform  load  was  applied  to  the  top  node. 
Fourteen  time  increments  were  chosen;  each  five  increments  represent  a 
tenfold  increase  in  time. 

The  finite  element  results  obtained  are  shown  in  Figures  4.6  and 
4.7,  plotted  along  with  the  theoretical  values  calculated  from  the  above 
formulae.  The  results  are  in  good  agreement.  More  accurate  results  can 
be  obtained  by  refining  the  part  of  the  mesh  near  the  load,  since  the 
excess  pore  pressure  gradient  is  large  in  this  area.  This  is  especially 
true  shortly  after  the  application  of  load.  Better  results  also  may  be 
obtained  by  taking  more  small  time  increments  at  the  start  of  the 
analysis.  Thus,  the  results  improve  at  greater  times. 
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Figure  4. 
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: Finite  Element  Mesh  for  One-Dimensional  Consolidation 


- 69  - 


Consolidation  Ratio,  Uz 


Figure  4.6 : Excess  Pore  Pressures  for  One-Dimensional  Consolidation 
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Figure  4.7:  Settlement  Versus  Time  for  One-Dimensional  Consolidation 


4 . 9 CONSOLIDATION  OF  A CLAY  LAYER  UNDER  A UNIFORM  STRIP  LOAD 

The  coupled  set  of  equations  which  governs  the  consolidation  process 
is  such  that  closed  form  solutions  have  only  been  obtained  for  problems 
with  very  simple  boundary  conditions.  The  case  of  a strip  load  on  an 
elastic  half  space  was  solved  by  Schiffmann,  Chen  and  Jordan  (1969).  A 
closed  form  solution  has  also  been  obtained  for  the  problem  of  a strip 
load  on  an  elastic  layer  witli  a smooth  base,  Gibson,  Schiffmann  and  Pu 
(1970).  However,  problems  in  the  field  do  not  often  approximate  the 
above  boundary  conditions,  and  thus,  these  solutions  can  rarely  be 
app 1 i ed . 

A more  useful  problem  is  that  of  consolidation  of  an  elastic  clay 
layer  under  a uniform  strip  load  with  a rough  rigid  impervious  base. 
While  no  closed  form  solutions  for  this  problem  exist,  a finite  element 
approach  presented  by  Christian,  Boehmer  and  Martin  (1972)  provides 
results  which  can  be  used  to  check  those  obtained  by  the  present 
program.  The  finite  element  mesh  used  is  shown  in  Figure  4.8,  and 
consists  of  seventy  Q8P4  elements  in  plane  strain.  The  QSP4  element  is 
an  isoparametric  eight  noded  quadrilateral  element  with  pore  pressure 
degrees  of  feedom  at  the  corner  nodes.  Water  is  allowed  to  drain  at  the 
surface,  while  the  base  is  impervious.  The  boundary  conditions  on  the 
edge  furthest  from  the  load  have  little  effect  on  the  solution.  As  a 
matter  of  preference  these  nodes  are  fixed  in  order  to  reduce  the  number 
of  degrees  of  freedom  in  the  problem.  Fourteen  time  increments  were 
chosen  in  the  same  manner  as  those  used  for  the  problem  of 
one-dimensional  consolidation. 
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The  settlement  results  are  plotted  against  non-dimensional  time  for 
the  load  center  and  the  load  edge  in  Figure  4.9.  These  are  compared  to 
one-dimensional  consolidation  settlement  results.  The  results  agree 
with  those  reported  by  Christian,  Boehmer  and  Martin  (1972). 

4.10  EFFECT  OF  ELEMENT  CHOICE  ON  SOLUTION  ACCURACY 

During  the  course  of  this  work  a number  of  different  finite  elements 
was  used,  depending  upon  the  application.  Hence,  some  effort  was 
expended  to  determine  the  most  superior  element.  All  elements  have  an 
isoparametric  formulation.  This  enabled  geometries  with  curved 

surfaces,  such  as  tunnels,  to  be  easily  modelled  by  higher  order 
elements.  The  following  elements  were  programmed: 

1.  Q4P4,  a four-noded  quadrilateral  element  with  a pore  pressure 
degree  of  freedom  at  each  node. 

2.  Q8P4,  an  eight-noded  quadrilateral  element  with  a pore  pressure 
degree  of  freedom  at  the  corner  nodes. 

3.  Q8P8,  an  eight-noded  quadrilateral  element  with  a pore  pressure 
degree  of  freedom  at  each  node. 

4.  L2P2,  a two-noded  link  element  with  a pore  pressure  degree  of 
freedom  at  each  node. 

5.  L3P2,  a three-noded  link  element  with  a pore  pressure  degree  of 
freedom  at  the  end  nodes. 

It  should  be  noted  that  the  L2P2  element  is  obtained  by  the 
degeneration  of  the  Q4P4  element  into  one  dimension.  Similarly,  if  the 

Q8P4  element  is  degenerated  into  one  dimension,  the  L3P2  element  is 
obtained.  This  study  focused  upon  the  quadrilateral  elements  since  the 
same  considerations  of  choice  may  then  be  applied  to  the  linear 
elements. 
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Figure  4.8:  Finite  Element  Mesh  for  Strip  Load  on  a Clay  Layer 
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Figure  4,9:  Settlement  versus  Time  for  Strip  Load  on  a Clay  Layer 


The  problem  of  one-dimensional  consolidation  was  approached  using 
both  the  Q8P8  and  Q8P4  elements.  The  same  numbers  of  elements  and  time 
steps  were  used.  It  is  found  that  small  oscillations  of  excess  pore 
water  pressure  exist  through  space  for  the  Q8P8  element  at  any  given 
time,  but  not  for  the  Q8P4  element,  (see  Figure  4.10.)  At  first,  this 
might  seem  surprising  since  there  are  more  pore  pressure  degrees  of 
freedom  for  the  case  using  the  Q8P8  element.  The  reason  for  this 
appears  to  be  the  matching  of  the  stress  and  pore  pressure  distributions 
across  the  elements.  Considering  the  Q8P8  element,  it  is  seen  that  the 
displacement  shape  functions  are  quadratic.  Thus,  the  stress  variation 
across  the  element  is  linear,  while  the  pore  pressure  variation  across 
the  element  is  quadratic.  This  mismatch  between  the  orders  of 
interpolation  for  stress  and  pore  pressure  is  probably  the  reason  for 
the  apparently  poor  showing  of  the  Q8P8  element.  The  Q8P4  element  has 
equal  orders  of  interpolation  for  stresses  and  pore  pressures.  For  this 
reason,  hybrid  elements  have  become  popular  for  coupled  problems,  for 
example  Osaimi  (1976). 

A comparison  between  the  Q8P4  and  Q4P4  elements  was  made  using  the 
strip  loading  problem  described  in  Section  4.9.  The  same  numbers  of 
elements  and  time  increments  were  used.  The  prediction  of  settlement 
for  the  two  elements  was  very  close  even  though  for  the  Q4P4  element 
there  were  many  fewer  displacement  degrees  of  freedom.  However,  the 
prediction  of  pore  pressures  using  the  Q4P4  element  produced  large 
oscillations  in  time  for  a given  point  in  space,  (see  Figure  4.11.) 
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Figure  4.10:  Oscillations  in  Space  for  the  Q8P8  Element 
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Sandhu,  Liu  and  Singh  (1977) 


studied  various  elements  for 


one-dimensional  problems.  They  found  serious  deficiences  for  elements 
such  as  the  Q8P8.  This  is  supported  by  the  one-  and  two-dimensional 
studies  presented  here,  and  thus,  the  Q8P4  element  was  chosen  as  being 
the  most  suited  for  the  tunnel  analyses. 

4.11  EFFECT  OF  INTEGRATION  RULE  ON  SOLUTION  ACCURACY 

During  the  marching  process  an  integration  rule,  defined  by  the 
value  of  (3,  was  used.  Values  of  B from  zero  to  one  may  be  chosen.  When 
B = 0.5  the  rule  corresponds  to  a trapezoidal  rule,  and  when  B = 0.67 
the  rule  corresponds  to  a parabolic  rule.  It  can  be  shown,  Booker  and 
Small  (1975),  that  for  B greater  than  or  equal  to  one  half,  the  solution 
is  unconditionally  stable. 

The  effect  of  B on  the  solution  for  the  excess  pore  pressure  during 
one-dimensional  consolidation  was  tested.  A mesh  of  five  Q8P8  elements 
was  used  and  the  results  are  shown  in  Figure  4.12.  It  can  be  seen  that 
for  values  of  B greater  than  0.5  the  solution  becomes  damped,  while  for 
values  of  B less  than  0.5,  oscillations  in  the  solution  become  large. 

As  a result  of  this  study,  B = 0.5  was  used  for  the  tunnel  analyses. 
It  was  also  found  that  accurate  solutions  could  be  obtained  while  using 
five  time  increments  to  cover  a tenfold  increase  in  time.  This 
procedure  was  also  used  in  the  tunnel  analyses. 
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Figure  4.11:  Oscillations  in  Time  for  the  Q4P4  Element 


= 0-67  B - 0-55  B - 0-50  3=0-45 
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Figure  4.12:  Effect  of  Integration  Rule  on  Solution  Accuracy 


Chapter  5 


FINITE  ELEMENT  FORMULATION  FOR  LARGE  DEFORMATIONS 
5.  1 INTRODUCTION 

For  structures  built  with  materials  such  as  steel  and  concrete,  it 
is  unusual  for  displacements  or  strains  to  become  large.  However,  in 
geotechnical  engineering  the  modulus  values  for  soil  may  be  small, 
thereby  producing  large  deformations.  For  example,  large  settlements 
are  often  observed  when  embankments  are  constructed  on  soft  clay  or 
peat.  In  this  report,  the  primary  application  is  to  that  of  shield 
tunneling;  in  this  case,  the  strains  may  range  from  small  to  large 
depending  upon  the  soil  conditions  and  construction  technique.  The  most 
likely  location  for  large  strains  is  in  the  area  of  the  tail  void  which 
is  created  as  the  shield  advances.  Here,  the  soil  is  temporarily 
unsupported  and  can  collapse  against  the  liner,  leading  to  large  strain 
effects.  Obviously,  this  is  most  likely  to  occur  where  the  soil  is  a 
soft  clay  or  a non-sel f supporting  medium  like  sand.  Grouting  behind  the 
liner  segments  is  designed  to  prevent  this  collapse,  but  if  the  soil  is 
weak  enough,  the  grouting  will  come  too  late. 

5.2  CHOOSING  A REFERENCE  FRAME 

To  include  the  effect  of  finite  deformations,  additional  terms  need 
to  be  added  to  the  stiffness  matrix.  In  order  to  compute  these  terms, 
it  is  necessary  to  distinguish  between  two  possible  reference  frames. 
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In  a Lagrangian  analysis,  the  coordinate  system  remains  fixed  in 


space.  This  method  leads  to  an  unsymmetric  stiffness  matrix,  Hibbitt, 
Marcal  and  Rice  (1970).  This  is  an  undesirable  result,  since  the 
solution  of  the  resulting  equations  requires  approx i mate  1 y twice  as  many 
operations  on  a digital  computer.  However,  if  the  stress-strain  law  is 
already  unsymmetric,  as  is  the  case  for  non-assoc i at i ve  plasticity,  then 
this  method  may  be  desirable.  Carter,  Booker  and  Davis  (1977). 

In  a Eulerian  analysis  the  coordinate  system  is  fixed  in  the 
material.  This  method  is  approximated  by  taking  small  load  increments 
and  updating  the  nodal  coordinates  after  each  increment.  This  method 
should  more  precisely  be  called  an  'updated  Lagrangian'  analysis. 
McMeeking  and  Rice  (1975)  used  this  approach  to  obtain  a symmetric 
stiffness  matrix,  and  this  is  the  method  used  here. 

5 . 3 CHOOSING  APPROPRIATE  MEASURES  OF  STRESS  AND  STRAIN 

When  the  finite  element  equations  are  derived  from  the  virtual  work 
equations,  care  must  be  taken  in  the  definition  of  stress.  The 
relationship  between  different  measures  of  stress  is  given  by  Yamada  and 
Wifi  (1977).  The  equilibrium  equations  are  most  conveniently  written  in 
terms  of  nominal  or  Lagrange  stress.  However,  a symmetric  stress 
measure  such  as  Cauchy  stress,  a,  or  Kirkoff  stress,  t,  is  more 
desirable.  Furthermore,  the  stress-strain  law  must  be  written  using  a 
stress  measure  that  is  frame  indifferent.  The  Jaumann  rate  of  Kirkoff 
stress  fulfills  this  requirement,  Malvern  (1969).  It  should  be  pointed 
out  that  other  rates  of  stress  also  exist.  The  Jaumann  rate  is  denoted 
by  a superior  star. 
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The  stress-strain  law  may  then  be  written  as 


* 

Tij  = Dijkl  eki 

where  the  deformation  rate  tensor  is  given  by 


e i j 


1 ' 
2 . 


dv  j 2>Vj  ' 

2)Xj  Sx  i - 


The  Cauchy  true  stress  increment  is 
* 

dtf  i j = Tjj  - ajj  ejj  + cr  j k Wki  + s j k k j 
where  the  skew  symmetric  spin  tensor  is  given  by 


dv  } 

dVj  ' 

<bx  j 

dx  i - 

(5.1) 


(5.2) 


(5.3) 


(5.4) 


The  decision  to  use  the  Jaumann  rate  of  Kirkoff  stress  in  the 
stress-strain  law  is  arbitary.  Thus,  the  Jaumann  rate  of  Cauchy  stress 
could  have  been  used  where 

V*  * 

(5.5) 


5.4  INITIAL  STRESS  STIFFNESS  MATRIX 

The  initial  stress  stiffness  matrix,  [Ks],  must  be  added  to  the 
conventional  stiffness  matrix  to  give  the  total  structure  stiffness 
matrix. 
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[K] 


l B ] T [ D ] [ B ) dV  + [ K s ] 


(5.6) 


V 

The  initial  stress  stiffness  matrix  was  derived  by  HcNeeking  and  Rice 
(1975),  and  is  given  by 


[K, 


( N ] T , ^ a i j [Nkl,j  - 2 [ B j ] T o'  ^ j (B^j)  } d V 


(5.7) 


where  [B,j]  is  expressed,  in  terms  of  the  shape  functions,  as 


[B}j]  — [ N 5 ] > j + [Nj],j 


(5.8) 


This  matrix  arises  from  two  sources.  The  variation  of  the 
strain-displacement  matrix  gives  one  term,  since  for  small  strains  this 
variation  is  zero.  This  is  often  referred  to  as  the  'initial  stability 
matrix',  Zienkiewicz  (1977).  The  other  term  comes  from  the  more  precise 
definition  of  stress.  The  initial  stability  contribution  is  dependent 
on  the  stress  level,  and  thus  an  eigenvalue  problem  may  be  solved.  Cook 
(1974) . 

|K0  + XK  s | = 0 (5.9) 

where  K0  is  the  stiffness  matrix,  excluding  the  initial  stability  terms, 
and  A is  the  eigenvalue  dependent  on  the  stress  level.  Thus,  this 
eigenvalue  problem  is  similar  to  a buckling  load  problem. 
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A FORTRAN  subroutine  which  calculates  this  initial  stress  stiffness 


matrix  for  the  case  of  plane  strain  is  given  in  Appendix  B.  The 
subroutine  computes  the  matrix  for  an  element  using  an  algorithm  derived 
by  Derbalian  (1978).  The  element  matrix  may  be  written  in  partitioned 
form  as 


[ k s I 


1 1 


% 

% 

r- 

r* 

v. 

l k i 2 ) ' 



> 

• l k 2 1 1 

l k 2 2 1 ' 

] J | dpdv 


where  the  submatrices  are  given  by 


1 

[kn]  = -(a 22  “ s-mHSvv)  ~ l Sum.) 

2 

1 

l k -j  2 ] = “ ""(Oil  + ® 2 2 ^ l ] ~ Ol  2 U Sun  ] + l Sw  1 ) 

2 

( k 2 1 ] = [ k i 2 ] T 
1 

t k 2 2 1 ~ ^2  2HS|jl(jl]  ~ O'  2 2 l Syy  1 

2 


and  where,  for  example 


mT 

. 


(5.10) 


(5.11) 


(5.  12) 


(5.13) 


(5.14) 


(5.15) 


The  computation  of  this  matrix  was  found  to  add  significantly  to  the 
cost  of  calculating  the  stiffness  matrix.  Hence,  the  program  option  to 
use  the  initial  stress  stiffness  matrix  should  be  used  with  caution. 
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5.5  UPDATING  STRESSES  AND  DISPLACEMENTS 


After  the  displacement  increment  has  been  calculated,  the  total 
stresses  must  be  updated  by  adding  the  incremental  stresses.  In  a large 
strain  analysis  this  stress  increment  must  first  be  modified  to  allow 
for  rotations.  This  is  achieved  using  the  following  equations  for  the 
plane  strain  case,  Derbalian  (1978): 

* 

dtfii  = Tii  ~ ~ 2Rai2  (5.16) 


* 

d (T  2 2 = *22  " p22e  + 2ROi2 


* 

dCJ  i 2 = T i 2 “ 1 2 ® + R(<Tl1  - CT  2 2 5 


* 

da 33  = T33  - a33e 


where  the  rate  of  dilation  is  given  as 


e = ei 1 + e22 


and  the  spin  is  given  by 


1 

'bV2 

dVi' 

2 

.dXi 

bx  2- 

(5.17) 


(5. 18) 


(5. 19) 


(5.20) 


(5.21) 


A FORTRAN  subroutine  which  calculates  this  true  stress  increment  for 
the  case  of  plane  strain  is  given  in  Appendix  B. 

At  the  end  of  each  increment  the  coordinates  of  the  nodal  points  are 
updated  by  adding  the  corresponding  displacement  increments.  This 
process  allows  the  reference  frame  to  move  with  the  material. 
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It  has  been  pointed  out,  Gibson,  England  and  Hussey  (1967),  that  for 


large  deformation  problems  the  velocity  of  the  pore  fluid  in  Darcy's  law 
needs  to  be  the  relative  velocity  of  the  fluid  to  the  solid.  By 
updating  the  coordinate  system  this  objective  is  achieved. 


5.6  RESIDUAL  LOAD  CORRECTION 

The  modulus  values  for  a Cam  clay  material  that  is  yielding  will  not 
remain  constant  throughout  the  load  increment.  Hence,  equilibrium  will 
not  exist  exactly  at  the  end  of  the  increment.  Stress  corrections  due 
to  the  rotation  of  elements,  which  were  discussed  in  the  previous 
section,  may  lead  to  further  deviations  from  the  equilibrium  state. 

To  correct  for  this  situation,  the  out-of-bal ance  forces  are 
calculated  at  the  end  of  each  increment,  and  then  applied  as  nodal  loads 
during  the  next  increment.  These  forces  may  be  obtained  from 


[ B 1 T { ct } dV  - {f  } 


V 


(5.22) 


where  { f } is  a vector  containing  the  total  nodal  loads. 

These  forces  should  be  small  compared  to  the  applied  loads  or 
support  reactions.  If  they  become  large,  the  solution  is  becoming 
unstable  and  smaller  load  increments  should  be  used. 
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5 . 7 VERIFICATION  OF  STRESS  CORRECTION  DUE  TO  ROTATIONS 


To  verify  the  stress  correction  due  to  finite  rotations,  one  element 
was  rotated  through  45°  in  the  x-y  plane.  The  element  was  prestressed 
in  the  x direction.  The  solution  was  compared  to  the  exact  solution 
obtained  using  the  Mohr's  circle  transformation.  Two  separate  computer 
runs  were  performed,  one  using  25  increments  of  rotation  and  the  other 
using  50  increments.  The  results,  shown  in  Table  5.1,  indicate  that  the 
program  yields  correct  results  which  come  closer  to  the  exact  solution 
as  more  increments  are  used. 


TABLE 

5.1 

Rotation  of 

Stresses 

0 XX 

CT  yy 

0 xy 

Before  rotation 

1 . 000 

0.000 

0.000 

Mohr  circle 

0.500 

0.  500 

0.  500 

50  increments 

0.506 

0.494 

0.512 

25  increments 

0.513 

0.487 

0.525 

5 . 8 LARGE  DEFORMATION  ONE  - 0 I MENS  I ONAL  CONSOLIDATION 

Finite  element  solutions  are  found  for  the  case  of  one-dimensional 
elastic  consolidation  with  large  deformations.  In  this  case  rotations 
do  not  occur  and  the  initial  stress  stiffness  matrix  was  not  computed. 
Hence,  the  only  correction  used  from  that  of  small  strain  theory  was  to 
update  the  reference  frame  after  each  increment. 
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The  mesh  used  consisted  of  five  L3P2  elements. 


Nineteen  time 


increments  were  used  with  the  full  load  applied  on  the  first  increment. 
The  final  theoretical  settlement  may  be  obtained  by  integrating  the 

strain 

do 

de  = — (5.23) 

E 

Let  the  final  stress  be  q,  and  the  initial  and  final  depths  of  the  clay 
layer  be  Ho  and  H respectively.  Thus, 

H q 

' dh  T do 

= — (5.24) 

J h J E 

H0  0 

Thus,  the  ratio  of  the  final  depth  to  the  initial  is  given  by 
H -q/E 

— = e (5.25) 

Ho 

Four  finite  element  runs  were  made  with  different  values  of  the 
ratio  q/E.  The  results  are  given  in  Table  5.2. 

The  percent  consolidation  is  plotted  against  the  time  factor  in 
Figure  5.1.  As  expected,  the  clay  layer  consolidates  more  quickly  at 
higher  strains.  This  can  be  attributed  to  the  progressively  shorter 
drainage  paths  that  are  obtained  during  the  consolidation  process  when 
the  final  settlements  are  large.  These  results  are  in  agreement  with 
the  finite  difference  results  of  Olson  and  Ladd  (1979),  and  the  finite 
element  results  of  Carter,  Small  and  Booker  (1977). 
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TABLE  5.2 


Final  Settlements  Computed  by  Different  Methods 


9 

% final 

% final 

% final 

% error 

- 

se  1 1 1 emen t 

sett  1 ement 

settl ement 

between 

E 

smal 1 strain 

large  strain 

finite 

FE  and  LS 

tlieo  ry 

theory 

e 1 ement 

0.00001 

0.001 

0.001 

0.001 

0.  1 

0.  1 

10.0 

9.5 

9.6 

1 . 1 

0.5 

50.0 

39.3 

40.9 

4.  1 

1.0 

100.0 

63.  2 

68.4 

8.2 

Mesri  and  Rokhsar  (1974)  found  that  small  strain  theory  gave 
reasonable  consolidation  results  up  to  strains  of  32  percent.  Because 
much  smaller  strains  are  expected  in  the  tunnel  analyses,  and  in  view  of 
the  good  results  obtained  using  just  the  updated  reference  frame,  the 
computation  of  the  initial  stress  stiffness  matrix  will  be  omitted  for 
subsequent  work  in  this  report. 

5.9  SUMMARY 

An  updated  Lagrangian  formulation  has  been  outlined  which  enables 
problems  involving  large  strains  to  be  solved  accurately  with  a 
symmetric  stiffness  matrix.  For  problems  with  intermediate  strains,  a 
less  expensive  option  is  available  in  which  only  the  nodal  coordinates 
are  updated  between  increments. 
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Figure  5.1:  Effect  of  Large  Strains  on  One-Dimensional  Consolidation 


Chapter  6 


CONSOLIDATION  ANALYSES  OF  TUNNELS 


6 . 1 INTRODUCTION 

Tunnels  are  constructed  in  many  different  shapes  and  sizes,  and  they 
may  pass  through  an  even  greater  variety  of  soil  profiles.  When 
analysing  a proposed  tunnel,  these  site  characteristics  would  be  used  in 
defining  the  finite  element  mesh  and  soil  properties.  For  the  purpose 
of  gaining  an  understanding  of  the  behavior  of  tunnels  due  to 
consolidation,  a typical  large  diameter  tunnel  in  a homogeneous  clay 
layer  was  analysed.  Parametric  studies  were  then  performed  to  determine 
the  effects  on  the  behavior  due  to: 


1.  Stiffness  of  the  liner 

2.  Size  of  the  tail  void 

3.  Soil  strength 

A.  Permeability  of  the  soil 

5.  Initial  pore  pressures  set  up  by  outward  shoving  at  the  tunnel 
face . 


This  is  by  no  means  a complete  list  of  the  factors  influencing  the 
behavior  of  a tunnel,  but  it  does  serve  to  cover  some  of  the  more 
important  ones.  The  last  item  is  directed  at  the  question  of  effects  of 
the  initial  pore  pressures  created  by  a pressurized  face  shield,  as 
noted  in  Chapter  1 1 . 
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6 . 2 PROBLEM  DEFINITION 


6.2.1  Geometry 

The  finite  element  mesh  used  is  shown  in  Figure  6.1.  The  soil  is 
modelled  using  66  eight-noded  isoparametric  quadrilateral  elements  with 
pore  pressure  degrees  of  freedom  at  the  corner  nodal  points  (Q8P4).  The 
tunnel  liner  is  modelled  using  8 eight-noded  isoparametric  quadrilateral 
el ements  (Q8) . 

The  tunnel  is  located  at  a depth  of  44  feet  (13.5  m)  to  the 
springline.  A rigid  base  is  assumed  at  a depth  of  88  feet  (27  m)  below 
the  surface.  The  outside  diameter  of  the  tunnel  is  26  feet  (8  m).  The 
finite  element  mesh  extends  horizontally  to  a distance  equal  to  11  radii 
from  the  center  of  the  tunnel.  This  will  ensure  that  the  boundary  has 
an  insignificant  effect  on  the  results. 

The  concrete  liner  is  assumed  to  be  20  inches  (0.5  m)  in  thickness, 
a value  typical  of  that  used  in  practice  for  this  size  tunnel. 

6.2.2  Soil  Conditions 

The  soil  is  assumed  to  be  a homogeneous  clay  layer  with  shear 
strength  increasing  with  depth  as  shown  in  Figure  6.2.  A desiccated 
layer  with  a uniform  shear  strength  is  assumed  near  the  surface.  The 
parameters  used  are  typical  of  Boston  Blue  Clay.  This  soil  was 
simulated  since  it  is  typical  of  soft  clays  in  urban  areas  and  since 
there  is  an  abundance  of  data  on  it. 
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Figure  6.1:  Finite  Element  Mesh  for  Tunnel  Analyses 


Two  shear  strength  profiles  below  the  desiccated  crust  were  used;  the 
first  corresponds  to  an  overconsolidation  ratio,  OCR,  of  one  and  the 
second  to  an  OCR  of  two.  The  coefficient  of  lateral  earth  pressure  at 
rest,  k0>  and  the  ratio  of  the  shear  strength  to  the  effective  vertical 
stress  are  dependent  on  the  OCR  and  are  given  in  Table  6.1,  based  upon 
results  reported  by  Ladd  et  al . (1977). 

Parameters  for  the  Cam  Clay  model  were  obtained  from  the  reference 
by  Carter,  Randolph  and  Wroth  (1979),  who  had  analysed  Boston  Blue  Clay 
for  an  embankment  problem  using  the  Cambridge  model.  These  parameters 
are  given  in  Table  6.2.  The  soil  is  assigned  a density  of  115  pcf  (18 
kN/m3 ) . 

When  air  pressure  is  not  used,  the  overstress  factor,  OFS,  may  be 
defined  as 

Pz 

OFS  = — (6.1) 

c u 

where  p2  is  the  total  vertical  insitu  pressure  at  a depth  equal  to  the 
depth  of  the  springline. 

The  OFS  has  been  used  in  practice  as  a parameter  for  predicting  the 
performance  of  a tunnel,  Broms  and  Bennermark  (1967).  When  0CR=1,  the 
value  of  OFS  is  7.3,  and  some  difficulty  may  be  expected  during 

construction.  This  could  include,  for  example,  a tendency  of  the  shield 
to  tilt.  When  0CR=2,  the  value  of  OFS  is  4.9,  in  which  case  the  tunnel 
should  be  constructed  without  unusual  difficulties. 
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Cu  (psf) 


Figure  6.2:  Shear  Strength  Versus  Depth 
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TABLE  6. 1 

OCR  Dependent  Parameters 

OCR  = 1 

OCR  = 2 

k0  0.55 

0.8 

cu/cr  v'  0 . 3 

0.45 

TABLE  6.2 

Boston  Blue  Clay  Parameters  for  Cam  Clay  Model 


K 

0.03 

A 

0.  15 

e c s 

1.74 

M 

1.20 

G/c  u 

74 

The  soil  was  assigned  an  isotropic  permeabilty  of  5 x 10'7  ft/sec 
(1.5  x 10'5  cm/sec).  The  actual  value  chosen  is  unimportant  since  it 
has  an  effect  of  just  scaling  the  time  dependent  results.  Results 
involving  actual  times,  rather  than  non-dimensional  times,  will  be  given 
in  Section  6.7. 

6.2.3  Boundary  and  Initial  Conditions 

The  water  table  is  assumed  to  be  at  the  surface.  Drainage  is 
allowed  to  occur  at  the  surface  and  at  the  base.  Thus,  the  base 
represents  a dense  sand  layer  or  permeable  bedrock.  The  tunnel  liner  is 
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assumed  to  be  impervious.  This  may  be  considered  true  provided  the 
permeability  of  the  liner  is  two  orders  of  magnitude  louer  than  the 
surrounding  soil,  Fitzpatrick  (1980).  It  is  difficult  to  achieve  this 
degree  of  waterproo f i ng  in  a liner,  but  it  is  possible  with  the  aid  of 
careful  grouting  techniques. 

Initially  there  is  a hydrostatic  pore  pressure  distribution.  The 
at-rest  stresses  in  the  soil  are  calculated  by  a simulated  gravity 
turn-on  with  the  appropriate  value  of  ko- 

6.2.  4 Liner  properties 

A concrete  tunnel  liner  is  used  with  a density  of  125  pcf  (20 
kN/m3).  The  weight  of  a concrete  liner  is  a significant  percentage  of 
the  weight  of  soil  removed  from  a tunnel  during  excavation,  and  hence  it 
must  be  included  in  the  analysis.  The  stiffness  is  calculated  using  a 
value  of  3000  ksi  (20700  kPa)  for  the  Young's  modulus  of  concrete.  The 
liner  is  assumed  to  consist  of  continuous  elements;  use  of  the  3000  ksi 
modulus  in  this  model  simulates  the  effect  of  a cast- i n-p 1 ace  liner.  If 
precast  segments  were  used,  their  ring  stiffness  would  be  less  than  that 
of  a continuous  liner  even  if  they  were  the  same  thickness  because  of 
the  effects  of  the  numerous  bolted  connections.  In  order  to  approximate 
the  reduction  in  stiffness  due  to  the  bolting  of  the  segments,  the 
modulus  is  reduced  by  dividing  by  a 'liner  stiffness  reduction  factor'. 
In  Section  6.4  four  different  liner  stiffness  reduction  factors  were 
investigated.  A reduction  factor  of  40  was  chosen  as  best  representing 
a typical  segmented  liner.  The  shortening  of  the  vertical  diameter  of 
the  tunnel  liner  was  used  as  a parameter  in  making  the  choice.  Values 
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of  this  parameter  from  the  finite  element  analyses  were  compared  with 
those  reported  by  Deere  et  al  (1969). 

The  permeability  of  a segmented  liner  is  typically  greater  than  that 
of  a cast- i n-p 1 ace  liner.  However,  a successful  grouting  program  should 
be  capable  of  making  up  the  difference. 

6.2.5  Construction  Performance 

The  construction  performance  for  the  cases  with  no  initial  excess 
pore  pressures  is  modelled  by  varying  the  tail  void  gap.  Thus,  a small 
gap  represents  good  control  of  the  shield.  Two  different  values  of  this 
tail  void  gap  were  used:  2.0  and  3.6  inches  (50  and  90  mm).  These  tail 
void  values  represent  the  distance  from  the  soil  as  it  leaves  the  tail 
of  the  shield  to  the  outside  of  the  liner  at  the  crown.  This  gap  is  not 
uniform.  For  example,  when  the  crown  gap  is  3.6  inches,  the  springline 
gap  is  2.0  inches  (50  mm)  and  the  invert  gap  is  0.8  inches  (20  mm).  The 
soil  closes  in  at  different  rates  around  the  liner.  Contact  with  the 
liner  is  assumed  to  take  place  at  the  same  time  for  each  point  around 
the  liner.  In  practice  this  may  be  the  time  when  the  grouting  is 
injected. 

In  the  case  analyzed  to  determine  the  effects  of  initial  excess  pore 
pressures,  it  is  assumed  that  a pressurized  face  shield  is  used  which 
has  caused  an  initial  outward  expansion  of  the  soil  away  from  the  shield 
face.  This  outward  expansion  actually  occurs  out  of  the  plane  of  the 
tunnel  cross-section  being  analyzed  and  thus  cannot  be  exactly  simulated 
with  the  present  plane  strain  approach.  In  lieu  of  an  exact  model,  the 
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initial  pore  pressures  for  this  analysis  are  created  by  an  expansion  of 
the  tunnel  opening  in  the  plane  of  the  finite  element  mesh,  as  described 
subsequently  in  this  chapter. 


6.2.6  Summary 

Three  different  parameters  were  varied  in  the  finite  element 
analyses.  Table  6.3  shows  the  parameters  used  for  the  six  analyses. 


TABLE  6.3 

Parameters  Used  in  Finite  Element  Tunnel  Analyses 


Run  Number 

Liner  Stiffness 
Reduction  Factor 

Crown  Tail 
Void  Gap 
(inches) 

OCR 

Initial 

Excess 

P.P. 

1 

40 

3.6 

1 

No 

2 

100 

2.0 

1 

No 

3 

40 

2.0 

1 

No 

4 

20 

2.0 

1 

No 

5 

10 

2.0 

1 

No 

6 

40 

2.0 

2 

No 

7 

40 

3.2 

1 

Yes 

6.3  BEHAVIOR  OF  A TYPICAL  TUNNEL 

The  results  shown  in  this  section  were  obtained  using  an  OCR  of  one 
for  the  soil  and  a tail  void  gap  of  3.6  inches  (90  mm)  from  the  liner  to 
the  crown  of  the  tunnel.  A liner  stiffness  reduction  factor  of  40  was 
used . 
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The  analysis  was  performed  using  63  time  increments  in  the  finite 


element  program.  Loads  may  be  applied  as  required  during  any  time 
increment,  however,  the  incremental  time  must  not  be  zero  in  any 
i ncrement . 

Increment  0 corresponds  to  the  initial  conditions  with  at-rest 
stresses  in  the  soil.  The  equivalent  nodal  loads  necessary  to  simulate 
the  excavation  of  the  tunnel  are  applied  during  increments  1 through  40. 
The  incremental  time  for  these  40  increments  is  kept  very  small  to 
ensure  undrained  behavior.  After  increment  15  the  soil  has  closed  the 
tail  void  and  the  liner  is  activated.  The  weight  of  the  liner  was 
included  during  increments  16  through  40.  The  construction  of  the 
tunnel  is  complete  after  increment  40.  The  large  number  of  increments 
used  are  required  because  of  the  non-linear  nature  of  the  Cam  clay 
model  . 

Consolidation  was  allowed  to  occur  during  increments  41  through  63. 
No  loads  were  applied  during  this  part  of  the  analysis.  The  time 
increments  are  increased  on  a logarithmic  scale,  so  that  each  time 
increment  is  approximately  1.6  times  larger  than  the  previous.  Hence, 
the  time  increment  increases  by  an  order  of  magnitude  every  five 
increments.  After  increment  63  the  excess  pore  pressures  had 
di ssi pated . 

In  discussing  the  results,  four  particular  stages  of  construction 
will  be  referred  to : 

1.  Initial  conditions  prior  to  any  construction. 
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2.  When  the  tail  void  has  closed,  just  before  the  liner  is 

acti vated . 

3.  After  all  excavation  forces  have  been  applied,  but  before 
consolidation.  This  represents  the  end  of  construction. 

A.  Final  conditions,  after  consolidation  is  completed. 


6.3.1  Settl ements 

Surface  settlement  profiles  before  and  after  consolidation  are  shown 
in  Figure  6.3.  The  maximum  surface  settlement  occured  above  the  tunnel 
crown  as  expected.  This  settlement  was  2.08  inches  (53  mm)  before 
consolidation.  This  increased  by  32  percent  to  2.74  inches  (70  mm) 
after  consolidation. 


Surface  settlement  profiles  have  been  found  in  practice  to 
approximate  the  error  function,  Schmidt  (1969).  For  field  data  the 
surface  settlement,  s,  at  a horizontal  distance,  x,  from  the  centerline 
of  the  tunnel  is  typically  approximated  by 


s = s m a x exp 


x2  'j 
2 i 2 - 


(6.2) 


where  i is  the  horizontal  distance  from  the  centerline  of  the  tunnel  to 
the  point  of  inflection  of  this  bell-shaped  curve.  Plotting  the 
logarithm  of  settlement  against  the  square  of  the  distance  from  the 
centerline  of  the  tunnel  the  data  should  fit  a straight  line,  and  i may 
be  obtained  from  the  slope  of  the  best-fit  line.  The  width  of  the 
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settlement  trough  may  be  expressed  in  non-dimensional  form  as  i/a,  where 
a is  the  tunnel  radius. 

The  error  function  is  plotted  in  Figure  6.4  for  the  surface  profiles 
both  before  and  after  consolidation.  The  data  may  be  seen  to  fit  a 
straight  line  very  well  for  both  cases.  The  width  of  the  settlement 
trough  increases  from  2.2  before  consolidation  to  3.1  after 

consolidation.  This  confirms  that  settlements  are  more  widespread  as  a 
result  of  consolidation,  Peck  (1969). 

The  volume  of  the  settlement  trough,  Vs,  may  be  conveniently 
obtained  from 

Vs  - 2.5  i smaK  (6.3) 

This  is  usually  expressed  as  a percentage  of  the  tunnel  opening 
V5 

v s = x 100%  (6.4) 

Tia2 

The  volume  of  the  settlement  trough  was  found  to  be  2.3  percent  prior  to 
consolidation.  This  increased  to  4.3  percent  after  consolidation, 
showing  the  effect  of  consolidation  to  be  significant.  Results  of  the 
error  function  analysis  will  be  compared  with  values  obtained  with  a 
different  soil  strength  and  field  data  in  Section  6.5. 
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Surface  Settlement  (inches)  Surface  Settlement  (inches) 


Figure  6.3:  Typical  Surface  Settlement  Profiles 


Figure  6.4:  Surface  Settlements  Plotted  as  an  Error  Function 
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Displaced  mesh  plots  before  and  after  consolidation  are  shown  in 
Figure  6.5.  The  displacements  are  scaled  by  a factor  of  30.  Before 
consolidation  there  is  a general  movement  of  the  soil  towards  the  tunnel 
opening,  and  the  springline  of  the  tunnel  settles.  As  consolidation 
takes  place,  there  is  a continued  movement  of  soil  above  and  below  the 
tunnel  touards  the  opening.  However,  the  tunnel  undergoes  continued 
ovalization  which  pushes  the  soil  outward  at  the  springline.  The 
springline  also  rises  to  its  original  level  due  to  the  dissipation  of 
negative  excess  pore  pressures  below  the  invert. 

For  one-dimensional  consolidation  the  degree  of  consolidation  may  be 
simply  given  as  the  compression  at  time  T divided  by  the  final 
compression.  However,  this  definition  is  not  satisfactory  when  both 
positive  and  negative  excess  pore  pressures  are  dissipating 
simultaneously  because  they  have  a cancelling  effect  on  each  other.  The 
excess  pore  pressure  dissipation  constant,  Ue»  is  introduced  where 


Ue 


dA 


(6.5) 


Thus,  the  degree  of  consolidation  may  be  defined  as 


Ue 

U = 1 (6.6) 

( U e) ma  x 


A non-dimensional  time  parameter,  T,  is  also  introduced 


kGt 

T = (6.7) 

H2y  w 
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where  k is  the  permeability  of  the  soil. 


G is  the  shear  modulus  at  the  depth  of  the  springline, 
t is  the  time  used  in  the  finite  element  program, 

H is  a measure  of  the  drainage  distance. 

In  this  case  the  drainage  distance  is  equal  to  the  depth  to  the 
spr i ng 1 i ne  . 

A plot  of  U versus  T is  shown  in  Figure  6.6.  The  shape  of  this 
curve  is  similar  to  curves  for  other  consolidation  problems.  Most  of 
the  consolidation  takes  place  between  T=0.005  and  T=0.1. 

The  maximum  surface  settlement  is  plotted  against  T in  Figure  6.7. 
The  surface  settlement  appears  delayed  when  compared  to  the  U versus  T 
curve.  This  is  due  to  the  rapid  dissipation  of  negative  excess  pore 
pressures  which  exist  below  the  invert  and  have  a short  drainage  path  to 
the  permeable  base.  These  negative  excess  pore  pressures  do  not 
contribute  to  the  settlement,  but  do  contribute  to  the  degree  of 
consolidation,  U.  On  the  other  hand,  the  positive  excess  pore  pressures 
near  the  springline  have  a long  drainage  path,  so  that  the  resulting 
settlement  is  delayed. 

6.3.2  Excess  Pore  Pressures 

The  predicted  distribution  of  excess  pore  pressures  is  shown  in 
Figure  6.8  for  two  stages  of  construction.  The  first  stage  is  just 
prior  to  tail  void  closure  and  the  second  stage  is  at  the  end  of 
construction,  before  consolidation.  The  maximum  absolute  values  of 
these  excess  pore  pressures  near  three  locations  are  shown  in  Table  6.4. 
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Figure  6.5:  Displaced  Mesh  Plots  for  a Typical  Tunnel 
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Figure  6.6:  Degree  of  Consolidation  Versus  Time  for  a Typical  Tunnel 


T 


Figure  6.7:  Maximum  Surface  Settlement  Versus  Time 
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Figure  6.8:  Excess  Pore  Pressure  Contours  for  a Typical  Tunnel 
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The  pore  pressure  contours  show  two  distinct  areas. 


The  first  is  a 


large  area  of  positive  excess  pore  pressures  extending  outward  from  the 
springline.  This  bulb  of  positive  excess  pore  pressures  grows  slightly 
from  the  first  stage  of  construction  to  the  end  of  construction.  The 
second  area  lies  beneath  the  invert  and  consists  of  very  high  negative 
excess  pore  pressures.  Although  the  extent  of  this  area  remains 
constant  as  construction  takes  place,  the  magnitude  of  these  pore 
pressures  increases  significantly.  The  pore  pressures  at  the  crown 
remain  almost  constant.  The  maximum  excess  pore  pressure  at  the 
springline  corresponds  to  a 20  percent  rise  in  the  hydrostatic  pore 
pressure,  while  the  excess  pore  pressure  at  the  invert  corresponds  to  a 
39  percent  drop  in  the  hydrostatic  pore  pressure. 


TABLE  6.4 

Excess  Pore 

Pressures  Near  Tunnel 

Crown 

Spr i ng 1 i ne 

Invert 

Before 

After 

Liner 

Liner 

Instal 1 at i on 
Instal 1 ati on 

10  psf 
- 1 0 psf 

509  psf 
549  psf 

-397  psf 
-1398  psf 

The  prediction  of  the  negative  excess  pore  pressures  in  the  crown 
and  invert  is  due  to  the  expansion  of  the  soil  into  the  tunnel  opening. 
Clough  and  Schmidt  (1977)  have  predicted  such  behavior  based  on  general 
principles.  The  very  high  values  of  negative  pore  pressure  predicted  in 
the  soil  below  the  invert  are  in  part  due  to  the  presence  of  the 
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underlying  rigid  base.  This  boundary  condition  creates  a situation 
where  the  soil  is  restrained  on  the  lower  side  while  it  tries  to  move 
into  the  opening  on  the  upper  side,  resulting  in  a stretching  effect. 
The  positive  pore  pressures  predicted  at  the  springline  are  generated 
because  this  area  is  subject  to  an  increase  in  both  shear  and  normal 
stress  as  the  tunnel  is  opened. 

6.3.3  Soil  Pressures  Acting  on  the  Liner 

The  in-situ  stress  distribution  at  the  location  of  a lining  prior  to 
the  excavation  of  a tunnel  is  elliptical.  The  pressure  at  the  imaginary 
crown  is  yz,  while  at  the  springline  it  is  koTz.  If  the  tunnel  is 
excavated  and  a perfectly  flexible  liner  installed,  the  liner  will 
deform  until  equilibrium  is  reached.  The  soil  pressure  distribution 
acting  on  a circular  liner  will  become  circular.  Direct  compressive 
stresses  will  exist  in  the  liner  and  there  will  be  no  bending  stresses. 
On  the  other  hand,  if  a rigid  liner  were  to  be  installed,  the  soil 
pressures  acting  on  the  liner  would  remain  unchanged  from  the  in-situ 
state . 

The  soil  pressures  acting  on  the  liner  before  and  after 
consolidation  are  shown  in  Figure  6.9.  The  in-situ  stress  state  is  also 
shown  for  comparison.  All  pressures  have  been  normalized  by  dividing  by 
the  vertical  in-situ  stress  at  the  springline.  It  can  be  seen  that  the 
elliptical  in-situ  distribution  becomes  more  circular  after  the  liner  is 
installed.  This  behavior  is  in  agreement  with  Peck's  (1969)  philosophy, 
outlined  above.  The  distribution  does  not  become  completely  circular 
because  the  tunnel  liner  is  not  perfectly  flexible.  Thus,  equilibrium 


is  maintained  by  the  bending  stresses  in  the  liner. 
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Figure  6.9:  Soil  Stresses  Acting  on  a Typical  Liner 


During  conso 1 i dat i on , the  excess  pore  pressures  near  the  springline 
decay,  reducing  the  horizontal  pressure  acting  on  the  liner.  The  decay 
of  negative  excess  pore  pressures  near  the  invert  causes  an  increase  in 
vertical  pressure  on  the  liner.  Thus,  the  distribution  of  soil  pressure 
acting  on  the  liner  becomes  more  elliptical  as  consolidation  proceeds. 

The  average  soil  pressure  acting  on  the  liner  is  found  to  increase 
by  3.1  percent  during  the  consolidation  period.  Much  larger  increases 
with  time  have  been  recorded  in  the  field.  Peck  (1969).  There,  the  soil 
pressure  acting  on  the  liner  is  often  calculated  from  strain  gauges 
which  are  attached  to  the  liner.  These  gauges  record  the  increase  in 
the  bending  stress  as  well  as  the  direct  stress.  It  is  possible  that 
the  average  stress  may  remain  constant  while  the  bending  stress 
i ncreases . 

The  elliptical  distribution  of  soil  pressure  on  the  liner  may  be 
approximated  by  two  terms  from  a Fourier  series: 

p = p0  + p2  cos  28  (6.8) 

Generally  the  average  pressure,  Po,  is  much  larger  than  the  other  term, 
p 2 • If  we  assume  that  the  liner  behaves  as  a thin  shell  with  the 
radius,  r,  much  greater  than  the  thickness,  t,  then  the  Po  term  produces 
a direct  compressive  stress,  given  by 

Por 

ad  = (6.9) 

t 
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No  bending  stress  is  produced  by  this  term. 


The  maximum  bending  stress 


produced  by  the  p2  term  may  be  obtained  from  the  Novozhilov  shell 
equations,  and  is  given  by 

2P2  r 2 

ffb  = (6.  10) 

t2 

The  direct  stress  for  the  p2  term  is  small  compared  to  the  bending 
stress . 

The  bending  component,  p2,  was  found  to  more  than  double  during 
consolidation,  as  shown  in  Figure  6.10.  This  phenomenum  is  one 

explanation  for  the  field  observations. 

The  horizontal  diameters  of  tunnel  liners  have  been  found  to 
increase  with  time  as  consolidation  takes  place.  Peck  (1969).  This 
effect  may  be  predicted  on  the  basis  of  the  excess  pore  pressure 
distribution  before  consolidation,  and  is  plotted  as  a percentage  of  the 
initial  diameter  versus  the  non-dimensional  time  parameter,  T,  in  Figure 
6.11.  The  increase  falls  within  the  range  of  those  recorded  in  the 
field. 

Increases  in  both  liner  pressures  and  liner  spreading  recorded  in 
the  field  tend  to  be  indefinite.  If  existing  tunnels  are  to  remain 
safe,  however,  these  increases  must  eventually  stop.  The  consolidation 
results  shown  in  Figures  6.10  and  6.11  indicate  that  these  increases  do 
indeed  die  out.  The  fact  that  this  has  not  been  recorded  in  the  field 
may  be  because  tunnels  are  not  monitored  long  enough  after  construction. 
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Normalized  Bending  Load, 


Figure  6.10:  Bending  Component  of  Pressure  on  a Liner  Versus  Time 
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Figure  6.11:  Typical  Liner  Spreading  as  a Function  of  Time 
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6.3.4  Soil  Stresses 


The  stress  paths  undergone  by  soil  elements  close  to  the  tunnel 
crown,  springline  and  invert  are  shown  in  Figure  6.12.  The  effective 
stress  paths  are  plotted  using  the  mean  effective  stress,  p',  and  the 
octahedral  shear  stress,  q.  The  p'  and  q axes  were  chosen  since  the  Cam 
clay  model  was  used  and  since  the  effect  of  02  is  important.  A basic 
reference  line  shown  on  this  diagram  is  the  critical  state  or  failure 
envelope,  defined  as 

q = Mp  (6.11) 

Figure  6.12  shows  that  the  soil  near  the  springline  moves  towards 
failure  during  excavation  and  away  from  failure  during  consolidation. 
The  soil  near  the  invert  behaves  in  an  opposite  manner,  that  is,  it 
moves  away  from  failure  during  excavation  and  towards  failure  during 
consolidation.  The  soil  near  the  crown  remains  almost  on  the  ko  line. 

Contour  plots  of  the  percentage  of  critical  state  reveal  the  areas 
moving  towards  failure.  The  situations  before  and  after  consolidation 
are  shown  in  Figure  6.13.  In  the  at-rest  state  prior  to  any  soil 
removal,  the  percentage  of  critical  state  is  53.  During  consolidation 
the  location  in  the  soil  that  is  nearest  to  failure  moves  from  near  the 
springline  to  near  the  invert. 

Relative  magnitudes  of  the  principal  stresses  and  their  directions 
are  shown  before  and  after  consolidation  in  Figure  6.14.  For  an  unlined 
tunnel  the  major  principal  stresses  would  be  aligned  tangent  to  the 
tunnel  for  a soil  element  next  to  the  tunnel  wall. 
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Figure  6.12:  Typical  Effective  Stress  Paths  for  Soil  Close  to  the  Liner 


Figure  6.13:  Percentage  of  Critical  State  Contours  for  a Typical  Tunnel 
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Figure  6,14;  Principal  Stresses  for  a Typical  Tunnel 
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In  this  case  the  stress  component  perpendicular  to  the  tunnel  wall  is 


carried  by  the  liner.  It  should  be  noted  that  the  direct  stress  in  a 
direction  parallel  to  the  axis  of  the  tunnel,  a2,  is  not  always  the 
intermediate  principal  stress,  a2.  For  example,  after  the  installation 
of  the  liner  the  stresses  at  a point  near  the  crown  were:  ax=1336psf, 
(Ty=1261psf,  a2=880psf  and  Txy=214psf.  Thus,  o2  is  actually  03.  This 
emphasizes  the  importance  of  including  ct2  in  the  finite  element 
calculations,  even  during  a plane  strain  analysis. 

6 . 4 EFFECT  OF  THE  LINER  STIFFNESS 

The  liner  used  in  the  finite  element  program  is  continuous,  while 
that  used  in  the  field  is  segmented.  The  bolted  joints  between  segments 
act  as  hinges  allowing  greater  movement.  In  order  to  simulate  this  more 
flexible  liner,  a reduction  factor  may  be  used  to  reduce  the  liner 
stiffness  calculated  for  one  which  is  continuous.  A suitable  liner 
stiffness  may  be  chosen  as  the  one  which  allows  liner  deformation 
similar  to  that  observed  in  the  field. 

The  effect  of  liner  stiffness  on  tunnel  behavior  is  investigated  by 
varying  the  elastic  modulus  for  the  concrete  liner.  Four  cases  were 
studied  with  modulus  reduction  factors  of  10,  20,  40  and  100.  The  liner 
became  more  flexible  as  the  reduction  factor  was  increased.  Other  basic 
parameters  for  the  finite  element  computer  runs  were  kept  constant.  The 
tail  void  gap  was  2.0  inches  (50  mm)  and  the  soil  was  normally 
consol i dated . 
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6.4.1  Settl ements 


Maximum  surface  settlements  are  shown  for  three  stages  of 
construction  in  Table  6.5.  The  maximum  surface  settlements  before  and 
after  consolidation  are  plotted  in  Figure  6.15. 


TABLE  6.5 

Variation  of  Settlements  with  Liner  Stiffness 


Reduction  Maximum  Surface  Settlement  (inches) 

Factor 


At  Tail  Void 

Before 

After 

C 1 osure 

Conso 1 i dat i on 

Consol i dat i on 

10 

0.90 

0.67 

0.38 

20 

0.90 

0 . S3 

0.73 

40 

0.90 

1.05 

1.23 

100 

0.90 

1 . 58 

2.24 

A flexible  liner  allows  more  movement  near  the  tunnel.  As  a result, 
a larger  zone  of  soil  undergoes  a change  of  stress.  When  this  soil 
yields  it  is  remolded.  The  larger  the  remolded  zone,  the  greater  is  the 
extent  of  excess  pore  pressure  buildup.  Greater  surface  settlements  are 
observed  after  consolidation  for  this  case. 

For  very  stiff  liners  the  settlement  actually  gets  reduced  from  the 
time  the  tail  void  is  closed  until  the  liner  is  installed.  This  result 
may  be  explained  by  examining  the  excavation  procedure,  equivalent  nodal 
loads  are  applied  to  the  nodes  on  the  tunnel  wall  to  simulate  the 
excavation.  These  loads  act  in  an  inward  direction.  After  the  crown 
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displacement  is  equal  to  the  tail  void  gap,  in  this  case  2.0  inches,  the 
liner  is  installed.  The  remaining  portion  of  the  equivalent  nodal  loads 
is  then  applied  to  the  outside  of  the  liner.  The  downward  load  at  the 
crown,  however,  is  less  than  the  upward  load  at  the  invert.  Thus,  there 
is  a net  upward  force  applied  to  the  liner.  When  the  liner  is  very 
stiff  it  acts  as  a rigid  body,  and  so  the  upward  force  is  translated 
into  an  upward  movement  at  the  surface.  During  the  consolidation  period 
the  settlement  is  further  reduced  because  of  the  dissipation  of  large 
areas  of  negative  excess  pore  pressures. 

In  order  to  identify  a reasonable  reduction  factor  for  modelling  a 
typical  tunnel,  the  shortening  of  the  vertical  diameter  was  compared  to 
those  reported  by  Deere  et  al . (1969).  The  percentage  change  in  this 
diameter  is  given  in  Table  6.6,  and  on  this  basis  a reduction  factor  of 
40  was  chosen  as  being  the  most  representative. 


TABLE 

6.6 

Effect  of  Liner  Stiffness  on 

Vertical  Diameter  Shortening 

Reduction  Factor 

Vertical  Diameter 

Shortening  (%) 

10 

0.4 

20 

0.6 

40 

0.9 

100 

1.4 
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6.4.2  Liner  Stresses 


The  p0  and  p2  components  of  pressure  before  and  after  consolidation 
are  shown  in  Table  6.7.  The  corresponding  stresses  are  plotted  in 
Figure  6.16.  The  direct  stress  increases  very  slightly  with 
consolidation,  while  the  bending  stress  approximately  doubles. 

Further,  while  the  direct  stress  is  only  very  slightly  reduced  by 
using  a more  flexible  liner,  the  bending  stress  experiences  a large 
reduct  ion. 


TABLE  6.7 

Effect  of  Liner  Stiffness  on  Liner  Pressures 

Pressures  as  a Percentage  of  the  Vertical  Insitu 
Stress  at  the  Springline 

Reduction  Before  Consolidation  After  Consolidation 
Factor 


Po 

P 2 

Po 

P 2 

10 

71.7 

9.4 

71 . 9 

18.2 

20 

71  . 1 

7.0 

71  . 6 

14.6 

40 

70.0 

5.5 

71.0 

11.5 

100 

67.3 

4.3 

69.3 

8.8 

6.4.3  Excess  Pore  Pressures  and  Soil  Stresses 

The  maximum  excess  pore  pressure  near  the  springline  increases  with 
an  increase  in  the  flexibility  of  the  liner.  These  pore  pressures  are 
shown  in  Table  6.8.  Their  increase  is  due  to  the  greater  movement 
allowed  by  a flexible  liner. 
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Figure  6.15:  Effect  of  Liner  Stiffness  on  Surface  Settlement 
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Liner  Stresses  (Ksi) 


0 10  20  40  100 


Liner  Stiffness  Reduction  Factor 


Figure  6.16:  Effect  of  Liner  Stiffness  on  Liner  Stresses 
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Excess  pore  pressure  contours  for  the  most  flexible  and  most  stiff  cases 


are  shown  in  Figure  6.17.  The  greater  extent  of  excess  pore  pressures 
is  clearly  shown  for  the  more  flexible  liner. 


TABLE 

CT» 

CO 

Effect  of  Liner  Stiffness  on 

Maximum 

Excess  Pore  Pressure 

Reduction  Factor 

Pore 

Pressure  (psf) 

10 

388 

20 

424 

40 

457 

100 

482 

The  corresponding  contour  plots  for  percentage  of  critical  state  are 
shown  in  Figure  6.18.  Again,  the  most  flexible  liner  produces  the 
greatest  yield  zone. 

6.4.4  Summary 

The  use  of  a stiff  liner  will  reduce  movements  close  to  the  liner. 
Thus,  the  buildup  of  excess  pore  pressures  will  be  minimized,  ultimately 
reducing  surface  settlements. 

The  use  of  a flexible  liner  will  greatly  reduce  the  bending  stresses 
in  the  1 iner . 
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Figure  6.17:  Effect  of  Liner  Stiffness  on  Excess  Pore  Pressures 
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Figure  6.18:  Effect  of  Liner  Stiffness  on  Percent  of  Critical  State 


6 . 5 EFFECT  OF  THE  TAIL  VO  I D SIZ E 


The  effect  of  construction  performance  is  modelled  by  varying  the 
size  of  the  tail  void  gap.  Two  values  of  this  parameter  were  used;  2.0 
inches  (50  mm)  and  3.6  inches  (90  mm).  The  2.0  inch  value  corresponds 
to  good  control  of  the  alignment  of  the  shield,  whereas  the  3.6  inch 
value  may  often  be  experienced  in  practice  when  less  than  ideal 
conditions  exist. 

6.5.1  Settl ements 

A comparison  of  maximum  surface  settlements  at  three  stages  of 
construction  is  given  in  Table  6.9.  The  settlements  are  considerably 
larger  when  construction  performance  is  poor.  The  increase  in 
settlement  due  to  consolidation  is  much  more  significant  when  larger 
initial  movements  are  allowed. 

For  the  Thunder  Bay  tunnel  discussed  in  Chapter  II,  settlements 
increases  due  to  consolidation  ranged  from  15  to  100  percent  of  the 
initial  settl ements . One  factor  contributing  to  the  larger  increase  was 
the  difficulty  experienced  in  properly  aligning  the  shield  which  led  to 
a larger  tail  void.  This  effect  is  observed  in  the  finite  element 
parametric  study. 

The  surface  settlement  profile  before  and  after  consolidation  is 
shown  in  Figure  6.19.  By  comparing  this  to  Figure  6.3  it  can  be  clearly 
seen  that  the  size  of  the  tail  void  has  a big  influence  on  the 
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TABLE  6.9 


Effect  of  Tail  Void  on  Maximum  Surface  Settlements 


Maximum  Surface  Settlements 
(inches) 

Tai 1 Void  Tai 1 Void 

2.0  inches  3.6  inches 


At  Tail  Void  Closure 
Before  Consolidation 
After  Consolidation 


0.90  1.91 
1.06  2.08 
1.24  2.74 


time-dependent  component  of  settlement.  The  corresponding  error 
function  plots  are  shown  in  Figure  6.20.  The  volume  of  the  trough  is 
reduced  from  4.3  percent  to  1.9  percent  as  the  tail  void  gap  is  reduced 
from  3.6  inches  to  2.0  inches.  The  width  of  the  trough  remains  almost 
unchanged . 

6.5.2  Soil  Pressures  Acting  on  the  Liner 

An  increase  in  soil  movement  prior  to  liner  installation  causes  a 
reduction  in  the  average  liner  pressure  both  before  and  after 
consolidation.  The  pressure  term  that  produces  bending  stresses  is 
initially  slightly  smaller,  but  after  consolidation  it  becomes  slightly 
larger  for  the  case  of  the  larger  tail  void.  The  liner  pressure  terms 
are  shown  in  Table  6.10. 
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TABLE  6. 10 


Effect  of  Tail  Void  on  Liner  Pressures 


Pressures  as  a Percentage  of  the  Vertical  Insitu 
Stress  at  the  Springline 

Tail  Void  Before  Consolidation  After  Consolidation 
(inches) 


Po 

P 2 

Po 

P z 

2.0 

70.0 

5.5 

71.0 

11.5 

3.  6 

66.4 

5.2 

68.5 

11.9 

6.5.3  Excess  Pore  Pressures  and  Soil  Stresses 

When  greater  movement  is  allowed  before  the  liner  comes  into  contact 
much  larger  excess  pore  pressures  are  mobilized,  (see  Figure  6.21). 
Figure  6.22  shows  the  larger  zone  of  yielding  clay  close  to  the 
springline  of  the  tunnel  for  the  case  of  the  larger  tail  void.  These 
percentage  of  critical  state  contours  were  plotted  at  the  time  when  the 
tail  void  had  just  closed. 

6.6  EFFECT  OF  SOIL  STRENGTH 

The  effect  of  soil  strength  on  the  tunnel  behavior  was  investigated 
by  comparing  the  response  for  two  different  values  of  the 
overconsolidation  ratio.  The  values  of  OCR  used  were  one  and  two. 
Other  basic  parameters  for  the  finite  element  computer  runs  were  kept 
constant.  The  liner  stiffness  reduction  factor  was  40  and  the  tail  void 
gap  was  2.0  inches  (50  mm). 
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Surface  Settlement  (inches)  2 Surface  Settlement  (inches) 


6.19:  Surface  Settlement  Profiles  for  a Tail  Void  of  2 inches 


Figure  6.20:  Error  Function  Plot  for  a Tail  Void  of  2 inches 
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Figure  6.21:  Effect  of  Tail  Void  Size  on  Excess  Pore  Pressures 
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Figure  6.22:  Effect  of  Tail  Void  Size  on  Remolded  Zone 
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6.6.1  Sett  1 ements 


Smaller  surface  settlements  would  be  expected  for  higher  strength 
soil.  Since  the  excess  pore  pressures  generated  are  also  considerably 
less,  smaller  settlements  due  to  consolidation  would  also  be  expected. 
The  values  of  maximum  settlement  are  shown  in  Table  6.11.  The  surface 
settlement  profiles  before  and  after  consolidation  are  shown  in  Figure 
6.23.  The  maximum  settlement  increases  slightly,  but  further  away  from 
the  tunnel  centerline  some  decrease  in  settlement  is  observed. 


TABLE  6.11 

Effect  of  OCR  on 

Maximum  Surface 

Set  1 1 ement 

Max i mum 

Surface  Settlements 
( i nches) 

OCR  = 1 

OCR  = 2 

At  Tail  Void  Cl osure 

0.90 

0.77 

Before  Consolidation 

1 . 06 

0.77 

After  Consolidation 

1.24 

0.83 

The  error  function  are  plotted  in  Figure  6.24  to  obtain  the 
settlement  trough  properties  shown  in  Table  6.12,  where  2 is  the  depth 
to  the  spr i ng 1 i ne . 
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Surface  Settlement  (inches)  Surface  Settlement  (inches) 


Figure  6.23:  Surface  Settlement  Profiles  for  OCR  = 2 


Figure  6.24:  Error  Function  Plot  for  OCR  = 2 
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TABLE  6.12 


Effect  of  OCR  on  Settlement  Trough 


OCR  = 1 OCR  = 2 


OFS 

7. 

3 

4. 

,9 

vs% 

1 . 

,9 

0. 

,8 

i /a 

3. 

, 0 

1 , 

, 8 

z/2a 

1 . 

, 7 

1 . 

.7 

The  volume  of  lost  ground  in  the  finite  element  studies  may  be 
compared  to  field  data  reported  by  Schmidt  (1969).  For  both  values  of 
OFS,  the  finite  element  studies  produced  ground  loss  close  to  the 
smallest  ground  loss  observed  in  the  field  for  the  corresponding  values 
of  OFS. 

This  indicates  good  agreement,  since  the  tail  void  gap  of  2.0  inches 
represents  good  control  of  the  shield,  and  since  ground  loss  due  to  face 
movement  is  not  considered  in  this  study. 

The  width  of  the  settlement  trough  for  the  case  of  0CR=1  is  large 
compared  with  data  reported  by  Schmidt  (1969).  The  one  field  case  with 
a wider  settlement  trough  was  for  a deeper  tunnel  in  Ottawa,  where  some 
consolidation  also  took  place.  This  confirms  that  settlements  may  be 
more  widespread  when  consolidation  takes  place.  Peck  (1969). 
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6.6.2  Soil  Stresses  Acting  on  Liner 


The  average  liner  pressure  did  not  change  during  the  consolidation 
period  and  was  found  to  be  63  percent  of  the  vertical  insitu  stress  at 
the  springline.  This  is  smaller  than  the  average  liner  pressure  of  71 
percent  observed  for  the  weaker  soil  case.  Thus,  the  soil  is  more 
capable  of  providing  support  when  it  has  a higher  strength. 

The  bending  stresses  are  reduced  by  a factor  of  almost  five  by 
increasing  the  strength  of  the  soil.  This  effect  may  be  expected 
because  of  the  axisymmetric  behavior  of  the  tunnel  in  the  stronger  soil. 

6.6.3  Excess  Pore  Pressures  and  Soil  Stresses 

Excess  pore  pressure  contours  are  shown  in  Figure  6.25.  The  excess 
pore  pressures  generated  were  considerably  smaller  than  those  for  0CR=1. 
The  maximum  positive  excess  pore  pressure  of  207  psf,  halfway  between 
the  springline  and  the  crown,  is  small  compared  to  437  psf  obtained  near 
the  springline  for  0CR=1.  The  shape  of  the  excess  pore  pressure 
contours  is  considerably  different  for  the  two  cases.  Positive  excess 
pore  pressures  occur  at  the  crown  for  0CR=2,  while  small  negative  excess 
pore  pressures  were  obtained  at  this  location  for  0CR=1. 

Using  0CR=2  the  insitu  stresses  were  calculated  with  ko=0.8.  Thus, 
the  initial  conditions  correspond  to  q/Mpr0 . 1 9 . Axisymmetric  conditions 
are  more  closely  approximated  around  the  tunnel  opening.  Percentage  of 
critical  state  contours  are  shown  at  tail  void  closure  and  before 
consolidation  in  Figure  6.26.  The  contours  show  considerable  changes 
from  those  previously  shown  in  Figure  6.13  for  0CR=1. 
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Figure  6.25:  Excess  Pore  Pressure  Contours  for  OCR 
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Figure  6.26;  Remolded  Zone  for  OCR 


Before  Consolidation 


Figure  6.27:  Principal  Stresses  for  OCR  = 2 
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The  principal  stresses  shown  in  Figure  6.27  reflect  this 
axisymmetric  behavior.  The  principal  stresses  have  rotated  much  more 
completely  at  the  crown  and  invert  than  those  obtained  for  the  0CR  = 1 
case . 

6 . 7 EFFECT  OF  SOIL  PERMEABILITY 

The  permeability  of  the  soil  has  a purely  non-dimensional  effect  on 
the  results.  All  results  that  are  a function  of  time  include  the 
permeability  parameter.  Thus,  the  value  of  the  permeability  does  not 
alter  any  of  the  stresses  or  displacements  reported.  This  is  based, 
however,  on  the  assumption  that  the  permeability  is  small  enough  that  an 
insignificant  dissipation  of  excess  pore  pressures  takes  place  before 
the  liner  is  installed  and  the  heading  sufficiently  advanced. 

Two  representative  soil  permeabilities  are  chosen  and  the 
time-dependent  effects  of  settlement  and  liner  spreading  are  plotted 
against  real  time  in  Figures  6.28  and  6.29  respectively. 

6 . 8 EFFECT  OF  INITIAL  EXCESS  PORE  PRESSURES 

As  noted,  this  case  was  analyzed  to  develop  an  understanding  of  what 
effects  in  time  dependent  response  might  be  created  by  initial  pore 
pressures  set  up  around  the  tunnel  due  to  the  use  of  a pressurized  face 
shield  machine.  This  situation  can  easily  develop  with  a slurry  or 
earth  pressure  balance  shield  as  illustrated  in  Chapter  II.  It  was  also 
pointed  out  earlier  that  the  present  analytical  approach  utilizing  plane 
strain  assumption  is  only  approximate  in  this  case  since  the  initial 
excess  pore  pressures  are  assumed  to  be  developed  by  in-plane  movement, 
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Figure  6.28:  Settlement  Versus  Time  for  Different  Soil  Permeabilities 
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Figure  6.29:  Liner  Spreading  Versus  Time  for  Different  Permeabilities 
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whereas  they  are  actually  developed  by  out-of-plane  movements  in  front 


of  the  shield.  Additional  studies  will  be  undertaken  at  a later  date 
using  the  axisymmetric  program  option  to  explore  these  effects  further. 
However,  the  present  approach  allows  for  a reasonable  assessment  of  the 
general  importance  of  these  conditions. 

The  process  of  creating  the  initial  excess  pore  pressures  involves 
expanding  the  circular  tunnel  under  undrained  plane  strain  conditions 
using  a trapezoidal  pressure  distribution  with  the  highest  pressure  at 
the  tunnel  invert  and  the  lowest  at  the  crown.  The  nonuniform  pressure 
distribution  is  used  to  reflect  the  fact  that  the  soil  is  assumed  to 
have  an  increasing  stiffness  with  depth,  and  the  stiffer  soil  is  more 
likely  to  allow  a higher  pressure  to  be  developed  during  the  shove  of 
the  shield.  The  pressures  used  are  chosen  to  produce  a set  of  outward 
peripheral  movements  around  the  tunnel  which  are  in  line  with  those 
observed  recently  at  the  site  of  an  earth  balance  shield  project  in  San 
Francisco  (Clough,  et  al.,  1982).  The  resulting  movements  are  depicted 
in  Figure  6.30.  As  would  be  the  case  in  the  field,  the  movements  are 
nonuniform  - 3.6  in.  (9.1  cm),  1.2  in.  (3.1  cm)  and  1.0  in.  (2.5  cm)  at 
the  crown,  springline  and  invert  respectively. 

The  procedures  for  modeling  the  remainder  of  the  tunnel  loading  are 
the  same  as  those  used  previously.  It  is  assumed  that  the  soil  is 
displaced  inwards  to  fill  the  tail  void,  whereupon  the  liner  is 
installed.  A tail  void  gap  of  3.2  in.  (8.1  cm)  is  used  at  the  crown. 
Note  that  the  soil  is  not  assumed  to  have  to  recover  the  entire  initial 
heave  before  displacing  inward.  The  reason  for  this  is  that  during  the 
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Figure  6.30:  Outward  Tunnel  Displacements  Designed  to  Induce  Initial 

Pore  Pressure 
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actual  out-of-plane  movements  in  front  of  the  shield  when  the  soil  at 


the  perimeter  of  the  shield  moves  away,  new  soil  takes  its  place  by 
moving  from  in  front  of  the  shield.  Thus,  soil  is  always  in  contact 
with  the  shield  and  can  only  move  inwards  the  amount  of  the  tail  void 
gap. 

The  excess  pore  pressures  developed  after  expansion  of  the  tunnel 
periphery  are  shown  in  Figure  6.31.  High  positive  values  developed  in 
the  areas  around  the  crown  and  the  invert.  At  the  springline  negative 
pore  pressures  are  predicted.  These  different  behaviors  reflect  the 
fact  that  in  the  crown  and  invert  areas  the  shear  stress  levels  are 
increased  over  their  original  at-rest  values,  while  at  the  springline 
area,  stress  levels  are  decreased.  Also  in  Figure  6.31  the  excess  pore 
pressures  after  closure  of  the  tail  void  are  depicted.  The  primary 
effect  of  this  movement  is  to  eliminate  the  negative  pore  pressures  at 
the  springline,  while  the  positive  pore  pressures  are  largely 
unaffected.  Importantly,  the  excess  pore  pressure  distribution  at  this 
stage  prior  to  consolidation  is  quite  different  than  that  for  simple 
tail  void  closure  (see  Figure  6.21). 

The  contours  of  percent  development  of  critical  state  for  the 
initially  outwardly  displaced  tunnel  after  tail  void  closure  (Figure 
6.32),  show  that  less  of  the  soil  strength  is  mobilized  as  compared  to 
the  case  of  simple  tail  void  closure  (see  Figure  6.22).  The  largest 
difference  occurs  at  the  springline  where  only  a maximum  of  60%  of 
critical  state  is  mobilized  in  Figure  6.32,  while  upwards  of  80-90%  is 
mobilized  in  Figure  6.22. 
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Figure  6.31:  Excess  Pore  Pressures  for  Case  of  Initial  Outward 

Movements 
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figure  6.32:  Contours  of  Percent  Critical  State  Mobilized 
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Surface  settlement  profiles  after  tail  void  closure  and  consolidation 


are  given  in  Figure  6.33.  These  profiles  are  drawn  assuming  that  the 
ground  surface  is  initially  level,  since  the  amount  of  initial  heave 
predicted  by  the  plain  strain  analysis  is  not  considered  representative 
of  the  actual  case.  Initially  the  settlement  trough  is  much  flatter 
than  those  observed  for  simple  tail  void  closure.  Subsequently,  after 
consolidation,  the  settlement  trough  takes  on  a more  conventional  form. 

The  effects  of  the  initial  heave  on  settlement  magnitudes  is 
generally  positive.  In  Table  6.13  the  maximum  movements  at  the  surface 
on  the  centerline  of  the  tunnel  are  tabulated  and  compared  to  those 
predicted  for  simple  tail  void  closure.  Even  assuming  no  initial  heave 
of  the  ground  surface  the  settlements  of  the  present  case  are  smaller  in 
every  instance  than  those  for  simple  tail  void  closure.  Since  the 
ground  surface  would  be  heaved  upwards  to  some  degree  by  out-of-plane 
movements,  the  resulting  settlements  would  be  even  smaller.  Thus,  for 
the  assumptions  made  in  these  analyses  and  levels  of  movement  assumed, 
it  appears  that  the  effects  of  the  initial  outward  displacement  caused 
by  a pressurized  face  shield  are  effective  in  reducing  movements.  The 
data  suggests  that  the  initial  heave  process  tends  to  "precompress"  the 
soil.  Actually  shear  stress  levels  in  the  springline  are  initially 
reduced,  and  shear  stress  increases  thereafter  cause  less  development  of 
excess  pore  pressures  during  subsequent  inward  movement.  Further  study 
of  this  problem  is  clearly  warranted  however,  to  determine  if  this 
finding  holds  for  a wide  range  of  conditions.  Very  large  initial 
outward  heaves  may  actually  be  detrimental  since  large  remolded  zones  of 
soil  could  be  created. 
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Surface  Settlement  (inches) 


Distance  from  Centerline,  x(feet) 


Figure  6.33:  Surface  Settlement  Profiles  for  Case  of  Initial 

Outward  Movements 
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TABLE:  6.  13 


Effect  of  Initial  Outward  Displacement  Caused  By  Pressurized  Face 

Sh i e 1 d 


Maximum  Surface  Settlement 

(inches) 

Simple  Tail  Void 

Initial 

Heave  Foil  owed* 

C 1 osure 

by  Tail 

Void  Closure  of 

2.0  Inch 

3.6  Inch 

3.2  Inch 

C 1 osure 

C 1 osure 

Closure 

At  Tail 

Void  Closure 

0.9 

1.91 

0.62 

Before 

Conso 1 i dat i on 

1.06 

2.08 

0.72 

After  Consolidation 

1.24 

2.74 

0.98 

*Note:  Assumes  Ground  Surface  Horizontal  After  Shove.  Any  Heave 

Which  Might  Occur  at  Surface  Would  Reduce  Settlements. 
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Chapter  7 


SUMMARY  AND  CONCLUSIONS 

Control  of  ground  movements  during  tunneling  is  one  of  the  major 
issues  related  to  construction  in  the  urban  environment.  New  shield 
tunneling  techniques  have  been  developed  to  address  this  problem  along 
with  improvements  in  conventional  methods.  However,  these  improvements 
have  not  been  paralleled  by  development  of  our  understanding  of  the 
reasons  for  the  behavior  of  the  ground  during  soil  tunneling.  In  the 
past,  movements  around  tunnels  have  largely  been  studied  by  field 
observations,  a particularly  effective  technique  where  the  data  can  be 
extrapolated  to  other  projects.  However,  these  procedures  cannot  be 
readily  used  to  predict  performance  for  new  tunneling  machines.  Also, 
since  measurements  are  often  discontinued  with  the  construction 
completion,  little  is  learned  about  long  term  effects,  which  in  some 
cases  can  be  very  important. 

This  report  describes  a research  program  devoted  to  the  study  of 
time-dependent  behavior  of  shallow  tunnels  in  clay  due  to  the  effects  of 
mobilization  and  dissipation  of  excess  pore  pressure.  To  the  knowledge 
of  the  authors  it  is  the  first  of  its  kind.  The  problem  is  investigated 
by : 

1.  Reviewing  available  field  data 
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2.  Developing  a finite  element  code  to  model  shield  tunneling 
effects  while  explicitly  taking  into  account  excess  pore  pressure 
effects  in  the  soil. 

3.  Performing  a series  of  parametric  studies  on  the  time-dependent 
behavior  of  tunnels. 

The  finite  element  program  which  is  developed  is  fully  described  and 
should  be  useful  in  future  studies  of  a range  of  advanced  tunneling 
prob 1 ems . 

As  expected,  the  review  of  existing  field  data  revealed  that  few 
cases  of  tunnel  behavior  have  been  monitored  long  enough  to  define  long 
term  time  effects.  Interestingly,  about  half  of  those  available  are 
very  recently  reported  and  concern  advanced  shield  tunnel  projects.  All 
of  the  data  available  shows  that  in  cohesive  soils  surface  settlements 
increase  with  time  following  tunnel  passage.  Most  authors  agree  that 
this  effect  can  largely  be  explained  as  due  to  dissipation  of  excess 
pore  pressures  set  up  in  the  soil  during  tunneling.  Also,  there  is 
strong  evidence  that  the  greater  the  degree  of  disturbance  of  ground 
during  tunneling,  the  greater  are  the  resulting  time-dependent 
movements.  In  addition  to  the  ground  movements,  liner  loads  and 

deformation  of  the  liner  continue  with  time.  Thus,  although  the  data 
are  somewhat  limited,  there  is  strong  evidence  for  significant  time 
effects,  particularly  in  the  case  of  the  newer  shield  types. 

In  order  to  study  the  time  dependent  behavior,  a finite  element  code 
was  developed  to  model  construction  of  a shield  tunnel  from  the  initial 
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condition  to  excavation  to  liner  installation,  and  finally  until  all 


excess  pore  pressures  have  dissipated.  In  particular,  the  finite 
element  code  was  designed  to 

1.  Simulate  the  construction  procedure. 

2.  Use  an  e 1 asto-p 1 ast i c soil  model. 

3.  Include  the  effect  of  consolidation. 

A.  Maintain  accuracy  for  large  strains. 

The  construction  procedure  was  simulated  in  plane  strain  since  the 
tail  void  is  the  most  significant  cause  of  settlements.  During 
excavation  the  soil  closes  in  on  the  opening 

by  an  amount  equal  to  the  size  of  the  tail  void  gap,  at  which  point  the 
liner  elements  are  activated.  The  liner  was  modelled  using  a ring  of  16 
eight-noded  isoparametric  elements  which  give  a good  response  in 
bending.  Undrained  conditions  were  assumed  during  the  construction 

sequence.  The  excess  pore  pressures  were  then  allowed  to  dissipate  and 

the  resulting  consolidation  behavior  was  observed. 

An  e 1 asto-p 1 ast i c soil  model  was  used  to  model  the  soil  yielding 
accurately  near  the  tail  void.  The  Cam  clay  soil  model  was  chosen  and 
developed  to  a form  suitable  for  coding  in  a finite  element  program. 
This  part  of  the  finite  element  code  was  verified  by  comparing  computed 
and  experimental  results  for  drained  triaxial  tests. 

\n  eight-noded  isoparametric  element  with  pore  pressure  degrees  of 
freedom  at  the  corner  nodes,  Q8P4,  was  used  to  model  the  soil. 
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Verification  of  the  consolidation  routines  was  performed  by  comparing 
computed  results  for  one-  and  two-dimensional  flow  and  consolidation 
problems  with  appropriate  closed  form  solutions. 

The  Q8P4  element  was  chosen  over  other  elements,  Q8P8  and  Q4P4, 
because  it  gave  more  accurate  results  for  the  above  problems.  The  time 
integration  parameter,  B,  was  set  equal  to  0.5  as  the  result  of  a 
stability  analysis  study. 

Large  strains  may  exist  in  the  soil  near  the  tail  void.  It  is 
important  that  accurate  solutions  be  obtained  in  this  area.  It  was 
found  that  sufficient  accuracy  could  be  obtained  by  updating  the  mesh 
geometry  after  each  time  increment  since  many  time  increments  are  used. 

Typical  tunnel  response  was  examined  using  the  finite  element  code 
for  a large  diameter  shallow  tunnel  in  Boston  Blue  Clay.  The  results 
obtained  are  in  agreement  with  reported  behavior  for  a tunnel  in  soft 
clay.  In  particul ar, 

1.  The  maximum  surface  settlement  was  found  to  increase 
significantly  during  consolidation.  The  width  of  the  settlement 
trough  also  increased. 

2.  The  bending  stresses  in  the  liner  were  found  to  increase 
significantly  with  time  although  the  average  stress  remained 
almost  constant. 

3.  Liner  spreading  was  observed  during  consolidation;  the  spreading 
stopped  once  the  excess  pore  pressures  had  dissipated. 

Parametric  studies  were  performed  to  determine  the  effects  on  the 
behavior  of  tunnels  due  to: 

1.  The  stiffness  of  the  liner. 
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2.  The  tail  void  size. 

3.  The  overconsolidation  ratio. 

4.  The  permeability  of  the  soil. 

5.  Initial  excess  pore  pressure  set  up  by  outward  pressures  at  the 
face  of  a pressurized  face  shield. 


The  following  is  a list  of  some  of  the  more  important  results  of  the 
parametric  studies: 


1.  A flexible  liner  allows  soil  movement,  which  leads  to  greater 
disturbance  and  thus  higher  excess  pore  pressures.  The  maximum 
surface  settlement  increases  more  during  consolidation  for  this 
case  than  for  a stiff  liner. 

2.  Bending  stresses  in  the  liner  are  reduced  as  the  liner  becomes 
more  flexible.  The  average  soil  pressure  acting  on  the  liner, 
however,  is  reduced  only  by  a small  amount. 

3.  Bending  stresses  in  the  liner  always  approximately  double  during 
consolidation  without  regard  to  liner  stiffness. 

4.  The  surface  settlement  increase  due  to  consolidation  becomes 
greater  as  the  tail  void  is  increased. 

5.  The  average  soil  pressure  acting  on*the  liner  decreases  slightly 
both  before  and  after  consolidation  with  an  increase  in  the  size 
of  the  tail  void. 

6.  The  behavior  of  a tunnel  located  in  overconsolidated  clay  is 
different  from  the  behavior  of  one  in  normally  consolidated  clay. 
In  normally  consolidated  clay,  the  maximum  excess  pore  pressures 
exist  near  the  springline  of  the  tunnel.  They  are  several  times 
greater  than  the  maximum  excess  pore  pressures  which  occur  near 
the  crown  of  a tunnel  for  an  overconsolidated  soil. 

7.  The  maximum  surface  settlement  remains  almost  unchanged  for  a 
soil  with  OCR  = 2. 

8.  The  average  soil  pressure  acting  on  the  liner  is  reduced  by  about 
12  percent  as  the  value  of  OCR  is  increased  from  one  to  two. 

9.  The  permeability  of  the  soil  is  contained  in  a non-dimensional 
time  parameter.  Thus  results  are  valid  for  a wide  range  of  soil 
permeabi 1 i ties. 

10.  Small  levels  of  initial  outward  displacements  around  the  tunnel 
periphery  due  to  pressure  at  the  face  produce  beneficial  effects 
in  reducing  time-dependent  movements  and  loadings. 
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Appendix  A 


PEPCO  USER'S  MANUAL 

A.  1 INTRODUCTION 

PEPCO  is  an  acronym  for  Program  for  E 1 asto-Pl ast i c consolidation. 
It  was  developed  at  Stanford  University  as  a research  tool  for  analysing 
the  consolidation  behavior  of  soil  using  the  finite  element  method.  The 
program  has  been  used  primarily  for  consolidation  analyses  of  tunnel 
behavior  in  clay.  Details  pertaining  to  the  theory  behind  the  program 
are  given  in  the  Ph.D.  report  entitled  'Finite  Element  Consolidation 
Analyses  of  Tunnel  Behavior  in  Clay'  by  Paul  R.  Johnston. 

The  program  is  coded  in  FORTRAN  IV  and  has  been  in  use  on  an  IBM 
3033  computer  at  Stanford's  Center  for  Information  Technology. 

A . 2 PROGRAM  CAPABILITIES 

The  capabilities  of  this  finite  element  program  are  controlled  by 
the  variety  of  element  and  material  types  in  its  respective  libraries. 
In  its  current  state  the  following  features  are  available: 

1.  Plane  strain  or  axisymmetric  analysis 

2.  One-dimensional  analysis 

3.  Elastic  analysis 

4.  E 1 asto-p 1 ast i c analysis  using  the  Cam  clay  soil  model 

5.  Fully  drained  or  undrained  analysis 
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6.  Fluid  flow  analysis 

7.  Consolidation  analysis 

8.  Large  strain  analysis 

9.  Initial  stress  generation 

10.  Installation  of  structural  members 

11.  Construction  of  embankments 

12.  Incremental  loading 

13.  Restart  capability 

14.  Expandable  element  library 

15.  Expandable  material  library 

16.  Output  control 

17.  Keyword  controlled  input 

It  is  relatively  easy  to  add  new  features  to  the  program  by  making 
additions  to  the  element  and  material  type  libraries.  For  example,  the 
program  may  be  converted  to  include  three-dimensional  problems  simply  by 
adding  a three-dimensional  element,  since  the  input,  solution  and  output 
routines  do  not  need  to  be  altered. 

A.  3 PROGRAM  OVERVIEW 

The  program  consists  of  a main  program  and  approximately  40 
subroutines,  see  Figure  A.1.  The  program  can  be  executed  using  a REGION 
of  512K  for  small  problems.  Double  precision  is  used  for  all  REAL 
variables.  The  array  storage  is  controlled  by  the  dimension  of  the 
vector  A in  the  main  program. 
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Figure  A. 1 : Program  Structure 


The  main  program  calls  the  subroutine  CONTRL. 


This  subroutine 


allocates  the  array  space  dynamically  by  mapping  the  smaller  arrays, 
both  INTEGER  and  REAL,  into  the  space  allocated  for  A.  CONTRL  then 
calls  the  other  main  subroutines  in  order. 

The  subroutines  RDCONT,  RDMESH,  RDINIT  and  RDLOAD  are  used  to  read 
the  input  data.  MS H GEN  is  a mesh  generation  routine  and  SVLOAO 
calculates  equivalent  nodal  loads  for  body  forces. 

The  subroutine  STIFF  assembles  the  structure  stiffness  matrix.  The 
subroutine  calls  the  appropriate  STIFnn  to  calculate  the  element 
stiffness  matrices.  The  element  stiffness  subroutines  call  other 
subroutines,  for  example,  to  calculate  the  shape  functions  or 
stress-strain  matrix. 

The  subroutine  SOLVE  solves  the  finite  element  equations  using  a 
skyline  solution  algorithm,  Mondkar  and  Powell  (1974).  The  bookkeeping 
for  the  stiffness  vector  is  performed  by  the  subroutine  DETNA. 

The  stresses  are  updated  using  the  subroutine  STRESS.  This 
subroutine  calls  the  appropriate  STREnn  for  each  element  to  calculate 
the  updated  element  stresses.  The  element  stress  recovery  routines  call 
other  routines,  for  example,  to  calculate  the  incremental  stress 
contribution  due  to  rotations. 

Results  are  printed  using  subroutine  OUTPUT.  They  may  also  be  saved 
on  a restart  file  using  WSTART.  A future  run  may  start  where  this  run 
left  off  using  RSTART . 
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ELEMENT  LI BRARY 


A.  4 

The  program  currently  has  eight  element  types.  New  elements  may  be 
added  to  the  program  by  including  two  new  subroutines  per  element. 

1.  STIFnn,  where  nn  is  the  element  type  number,  is  a subroutine  that 
calculates  an  element  stiffness  matrix. 

2.  STREnn  is  a subroutine  that  computes  the  incremental  element 
stresses  and  updates  the  total  stresses. 

The  subroutine  CONTRL  must  be  modified  to  include  the  following 
information  regarding  the  new  element  type: 

1.  Number  of  coordinates  per  node. 

2.  Maximum  number  of  degrees  of  freedom  per  node. 

3.  Number  of  integration  points  per  element. 

4.  Number  of  nodal  points  per  element. 

A figure  is  given  for  each  element  type  showing  the  degrees  of 
freedom  at  each  node,  the  node  numbering  order  for  the  element 
connectivity  list,  and  the  local  and  global  coordinate  systems. 

The  element  types  are: 


1.  A four-noded  isoparametric  quadrilateral  element,  Q4.  The 
stiffness  matrix  is  evaluated  at  four  integration  points.  The 
element  uses  the  following  shape  functions: 


( 1-f ) (1-7))/4 

(A.  1) 

( 1 + f ) ( 1 ~7)  ) /4 

(A. 2) 

( 1 + * ) ( 1+7))/4 

(A.  3) 

( 1-f)  ( 1+7))/4 

(A. 4) 
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Figure  A. 2:  Q4  Element  Type  1 


y 


u,  v,  p 


Figure  A. 3:  Q4P4  Element  Type  2 
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2.  A four-noded  isoparametric  quadrilateral  element  with  a pore 
pressure  degree  of  freedom  at  each  node,  Q4P4.  The  shape 
functions  for  both  pore  pressure  and  displacement  degrees  of 
freedom  are  identical  to  those  of  the  previous  element.  Four 
integration  points  are  used. 

3.  An  eight-noded  isoparametric  quadrilateral  element,  Q8.  The 
stiffness  matrix  is  evaluated  at  nine  integration  points,  and  the 
following  shape  functions  are  used: 


fl 

- 

-(  1-f  ) ( 1-7))  ( 1 + f + 7))/4 

(A.  5) 

*2 

= 

- ( 1 + ?)  ( 1-7))  ( 1-f  + 7))/4 

(A.  6) 

f3 

= 

-(  1 + f ) ( 1+7))  ( 1-f-7))/4 

(A.  7) 

u 

= 

-(  1-f)(1+7j)  ( 1 + f-7))/4 

(A. 8) 

u 

- 

( 1-f)  ( 1-7))  ( 1 + f)/2 

(A.  9) 

u 

= 

( 1 + f ) ( 1-7))  ( 1+7))/2 

(A. 10) 

*7 

- 

( 1 + f ) ( 1+7))(1-f)/2 

(A.  1 1) 

u 

= 

( 1 — f ) ( 1+7))  ( 1 7)  ) / 2 

(A.  12) 

4.  An  eight-noded  isoparametric  quadrilateral  element  with  pore 
pressure  degrees  of  freedom  at  the  corner  nodes,  Q8P4.  The 
displacement  shape  functions  are  the  same  as  those  used  in 
element  3,  while  the  pore  pressure  shape  functions  are  the  same 
as  those  used  in  element  1.  Four  integration  points  are  used. 
This  element  is  recommended  for  two-dimensional  consolidation 
anal yses. 
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Figure  A. 5:  Q8P4  Element  Type  4 
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5.  An  eight-noded  isoparametric  quadrilateral  element  with  pore 


pressure  degrees  of  freedom  at  all  nodes,  Q8P8.  The  shape 
functions  used  for  both  the  displacement  and  pore  pressure 
degrees  of  freedom  are  the  same,  and  they  are  identical  to  those 
used  in  element  type  3.  Four  integration  points  are  used. 

6.  A two-dimensional  beam  element,  BM2D . This  element  is  useful  for 
a plane  frame  type  of  analysis,  or  as  a structural  member  in  a 
plane  strain  geotechnical  problem.  The  stiffness  matrix  for  this 
elastic  beam  is  written  in  closed  form.  The  displaced  shape  is 
given  by 


A beam  element  must  not  share  a node  with  an  element 
possessing  a pore  pressure  degree  of  freedom. 

7.  A three-noded  isoparametric  link  element  with  pore  pressure 
degrees  of  freedom  at  the  end  nodes,  L3P2.  The  shape  functions 
for  the  displacement  degrees  of  freedom  are 


7)  = a0  + + a 2f2  + a3f3 


(A. 13) 


U = -f(1-f)/2 


(A. 14) 


f2  = f(1+*)/2 


(A. 15) 


f3  = W2 


(A. 16) 
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Figure  A. 6:  Q8P8  Element  Type  5 


Figure  A. 7 : BM2D  Element  Type  6 
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The  shape  functions  for  the  pore  pressure  degrees  of  freedom  are 


f i = ( W )/2  (A. 17) 

f 2 = ( 1+f )/2  (A. 18) 

The  element  has  two  integration  points.  It  has  limited  uses  but 
is  extremely  efficient  for  one-dimensional  consolidation 
prob 1 ems . 

8.  A tuo-noded  link  element  with  pore  pressure  degrees  of  freedom  at 
both  nodes,  L2P2.  This  element  was  developed  for  problems  of 
combined  axial  symmetry  and  plane  strain;  for  example,  the 
expansion  of  a cylindrical  cavity.  Two  integration  points  are 
used  and  the  displaced  shape  is  given  by 

u=cir+C2/r  (A. 19) 
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Figure  A. 8:  L3P2  Element  Type  7 
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Figure  A, 9:  L2P2  Element  Type  8 
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A . 5 MATERIAL  LIRRARY 


Materials  may  be  defined  as  elastic  or  e 1 asto-p 1 ast i c . The 

e 1 asto-p 1 ast i c material  is  defined  using  the  Cam  clay  model  developed  at 
Cambridge  University.  However,  because  it  is  necessary  to  make  a 
distinction  between  one-  and  two-dimensional  materials,  and  because  the 
beam  element  has  its  own  material  type,  five  material  types  exist.  As 
many  material  numbers  may  be  used  as  needed  to  descibe  the  problem;  for 
example,  five  kinds  of  Cam  clay  could  be  used  along  with  two  elastic 
materials  and  nine  different  beam  sections  for  a total  of  16  material 
numbers . 

New  material  types  may  be  added  to  the  program  by  expanding 
subroutines  DMAT1D,  D MAT  2D  and  UPDMAT.  The  subroutine  CONTRL  must  be 
modified  to  include  the  following  information  regarding  the  new  material 
type  : 

1.  Number  of  material  properties  needed  to  define  the  material. 

2.  Number  of  stress  and  strain  variables  to  be  stored  at  each 
integration  point. 

The  material  types  are: 

1.  Linear  elastic  material  for  two-dimensional  elements,  EL2D.  The 
material  properties  are  given  in  Table  A.1.  The  first  seven 
parameters  may  be  set  to  zero  when  elements  that  do  not  possess 
pore  pressure  degrees  of  freedom  are  being  used. 

The  input  for  the  first  three  parameters  may  be  explained 
using  an  example:  if  the  hydrostatic  head  increases  in  the 

negative  y direction  then  igx=0.0,  igy=-1.0,  and  ig2=0.0. 
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TABLE  A. 1 


Parameters  for  Material  Type  1,  EL2D 


Variable 


Expl anat i on 


PRC  1 , NMAT) 

1 gx 

X 

component 

of 

the 

hydrostat i c 

head 

PR (2, NMAT) 

1 sy 

y 

component 

of 

the 

hydrostat  i c 

head 

PRO, NMAT) 

igz 

z 

component 

of 

the 

hydrostatic 

head 

PR ( 4, NMAT ) 

the 

density 

of 

water 

PR ( 5 , NMAT ) kx  the  permeability  in  the  x direction 

PR(6,NMAT)  ky  the  permeability  in  the  y direction 

PRO, NMAT)  k 2 the  permeability  in  the  z direction 

PR(8,NMAT)  E Young's  modulus 

PR ( 9 , NMAT ) v Poisson's  ratio 


The  output  of  stress  and  strain  variables  at  each  integration 
point  is  as  foil ows : 

1.  e x,  the  strain  in  the  x direction 

2.  € y,  the  strain  in  the  y direction 

3.  e z,  the  strain  in  the  z direction 

4.  yXy»  the  shear  strain 

5.  ox',  the  effective  stress  in  the  x direction 

6.  CTy',  the  effective  stress  in  the  y direction 

7.  ex  z ' > the  effective  stress  in  the  z direction 

8.  t x y > the  shear  stress 

9.  u,  the  pore  pressure 

2.  Cam  clay  material  for  two-dimensional  elements,  CC2D.  This 
material  may  be  specified  with  a shear  strength  that  varies  with 
depth.  The  first  seven  properties  are  identical  to  those  for 
material  type  1.  The  remainder  of  the  parameters  are  given  in 

Tabl e A. 2. 
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TABLE  A. 2 


Parameters  for  Material  Type  2,  CC2D 


Variable  Explanation 


PR ( 8, NMAT ) 
PR  C 9 , NMAT ) 
PRC  10, NMAT) 
PRC  1 1 , NMAT) 
PRC  12, NMAT) 

PRC  13, NMAT) 
PRC  14, NMAT) 


K slope  of  the  rebound  isotropic 
consolidation  curve 
X slope  of  the  virgin  isotropic 
consolidation  curve 

ecs  the  voids  ratio  for  unit  p'  on  the 
critical  state  line 

M the  slope  of  the  failure  line  in  p'-q 
space 

m where  the  shear  modulus  G = msu 
if  the  value  is  less  than  one  half  the 
value  is  assumed  to  be  Poisson's  ratio  v 
s0  where  the  shear  strength  is  given 
by  su  = s0  + Sid  where  d is  the  depth 
Si  given  above 


Stresses  are  output  in  the  same  order  as  those  for  the 
previous  material  In  addition,  the  following  variables  are  output 
at  each  integration  point: 


10. 

Pc' 

, the  effective 

preconsolidation  stress 

1 1 . 

P'» 

the  effective 

mean  stress 

12. 

q» 

the  octahedral 

shear  stress 

13. 

e , 

the  voids  ratio 

3.  Material  properties  and  cross-sectional  geometry  for  a beam 
element,  BM2D.  The  first  seven  fields  are  left  blank.  The 
parameters  are  given  in  Table  A. 3.  The  six  member  end  actions 
are  printed  for  each  element. 
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TABLE  A. 3 


Parameters  for  Material  Type  3,  BM2D 
Variable  Explanation 


PR ( 8 , NMAT ) 
PR ( 9 , NMAT ) 
PRC  10, NMAT) 


E Young's  modulus 
A the  cross-sectional  area 
I the  second  moment  of  area 


4.  Linear  elastic  material  for  one-dimensional  elements,  EL1D.  The 
input  is  identical  to  that  for  material  type  1.  However,  since 
there  is  no  shear,  only  seven  of  the  variables  are  printed  at 
each  integration  point. 

5.  Cam  clay  material  for  one-dimensional,  elements,  CC1D.  The  input 
is  identical  to  that  for  material  type  2.  However,  since  yxy  and 
Txy  are  zero,  only  eleven  of  the  variables  are  printed  at  each 
integration  point. 


A. 6 INPUT 

The  input  for  PEPCO  is  separated  into  four  blocks,  as  shown  in 
Figure  A. 10,  and  read  from  unit  5.  Each  block  is  terminated  by  an 
appropriate  END  card. 
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Load  Block 
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Within  each  block  there  may  be  a number  of  groups  of  cards. 


Each 


group  of  cards  has  a descriptive  keyword  on  its  first  card.  For 
example,  the  first  card  in  the  material  properties  group  has  the  keyword 
MATERIAL  in  the  first  eight  columns,  and  the  cards  following  the  keyword 
card  contain  the  material  properties.  The  groups  of  cards  may  be 
arranged  in  any  order,  within  any  given  block.  Only  those  groups 
pertinent  to  the  problem  being  analysed  need  be  included.  For  example, 
if  one  does  not  wish  to  use  a restart  file,  the  RESTART  group  need  not 
be  used. 

A. 6.1  Control  Block 

This  block  contains  information  regarding  the  type  of  analysis,  the 
size  of  the  problem  and  the  number  of  load  increments  desired.  The 
information  provided  in  this  block  is  used  to  allocate  the  array  space. 

a)  Title  group 

This  group  of  cards  is  used  to  input  a descriptive  title  for  the 
problem.  The  group  may  be  repeated  if  more  than  one  title  line  is 
desi red . 


Title 

Card 

Col umns 

Var iabl e 

Format 

Expl anation 

1 

1 - 8 

WORD 

A8 

Enter  the  word  TITLE 

2 

1 - 80 

TITLE 

10A£ 

Enter  a title  for  the  problem 

b)  Size  group 
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This  group  of  cards  is  used  for  the  purpose  of  internally 


dimensioning  the  arrays. 


Problem  size 

Card  Columns  Variable  Format  Explanation 


1 1 - 

8 

WORD 

A8 

Enter 

the 

word  SIZE 

2 1 - 

5 

NN 

I 5 

Number 

of 

nodes 

6 - 

10 

NE 

15 

Number 

of 

e 1 ements 

11  - 

15 

NM 

15 

Number 

of 

materials 

16  - 

20 

NBC 

I 5 

Number 

of 

boundary  conditions 

c)  Large  displacement  option 

This  card  is  used  if  a large  displacement  analysis  is  to  be 
performed.  When  this  option  is  invoked,  the  initial  stress  stiffness 
matrix  is  added,  and  rotations  are  considered  when  updating  the 
stresses.  The  nodal  coordinates  are  updated  regardless  of  whether  this 
option  is  used. 


Large 

d i sp 1 acements 

Card 

Co  1 umns 

Var i ab 1 e 

Format 

Expl anation 

1 

1 - 8 

WORD 

A8 

Enter  the  words  LARGE  DI 

d)  Residual  load  correction 


This  option  may  be  used  with  or  without  the  large  displacement 
option.  However,  if  the  large  displacement  option  is  being  used,  this 
option  is  recommended  in  order  to  maintain  equilibrium.  By  using  this 
option,  the  out-of-balance  forces  are  calculated  at  the  end  of  each 
increment  and  applied  as  nodal  loads  during  the  next  increment. 


I 

e)  Print  group 

This  group  may  be  used  to  suppress  output  or  define  the  number  of 
load  increments  between  printings.  This  group  may  be  placed  in  any  of 
the  four  blocks.  Hence,  it  is  possible  to  turn  printing  on  and  off  as 
desired.  Users  are  warned  that  if  the  default  printing  is  used,  large 
amounts  of  output  are  produced.  Output  is  written  on  unit  6. 


Print 

Card 

Col umns 

Variable 

Format 

Expl anation 

1 

1 - 8 

WORD 

A8 

Enter  the  word  PRINT 

2 

1 - 5 

I PRT 

15 

Results  will  be  printed  every 
IPRT  increments. 

IPRT=1  is  the  default. 

Set  I PRT  = 0 to  suppress  output. 

f)  Restart  group 


Residual  load  correction 

Card  Columns  Variable  Format  Explanation 

1 1-8  WORD  A8  Enter  the  word  RESIDUAL 
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This  group  is  used  to  save  the  results  of  a run  on  a restart  file. 
The  option  may  then  be  used  subsequently  to  read  this  file  and  continue 
an  analysis.  If  the  read  option  in  this  group  is  used,  then  the  mesh 
block  and  initial  stress  block  should  be  omitted,  since  this  information 
was  provided  in  a previous  run  and  saved  on  the  restart  file.  This 
option  is  particularly  useful  in  consolidation  analyses  when  the  time 
taken  for  the  pore  pressures  to  dissipate  is  unknown.  It  is  also  useful 
in  e 1 asto-p 1 ast i c runs  in  order  to  break  up  a large  run  into  a number  of 
smaller  runs  and  reduce  the  chances  of  wasting  a great  deal  of  money  on 
one  bad  run.  Restart  files  are  written  without  format  statements  to 
unit  8,  and  they  are  similarly  read  from  unit  9. 


Restart 

Card 

Co  1 umns 

Variable 

Format 

Exp  1 anat i on 

1 

1 - 8 

WORD 

A8 

Enter  the  word  RESTART 

2 

1 - 5 

IRST 

15 

Increment  IRST  will  be  read 
the  restart  file. 

Set  IRST  =— 1 if  this  option 
requ i red . 

f rom 
is  not 

5-10 

I WST 

15 

Results  will  be  saved  on  a 
file  beginning  at  increment 
Set  I WST  =- 1 if  this  option 
required. 

restart 
I WST . 
is  not 

g)  Number  of  increments  group 


This  group  should  be  used  if  more  than  one  increment  is  desired. 
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I ncrements 


Card 

Co  1 umns 

Variable 

Format 

Exp  1 ana t i on 

1 

1 - 

8 

WORD 

A8 

Enter  the 

word 

INCREflEN 

2 

1 - 

5 

NFL  I 

15 

Number  of 

the 

first  1 oad 

i ncremen  t 

6 - 

10 

NLLI 

15 

Number  of 

the 

last  load 

i ncrement 

h)  Pore  pressure  degree  of  freedom  group 

This  group  must  be  used  to  define  the  time  integration  parameter 
when  pore  pressure  degrees  of  freedom  are  used  for  consolidation  or 
fluid  flow  prob 1 ems . 


Pore 

pressure 

Card 

Co  1 umns 

Vari abl e 

Format 

Explanation 

1 

1 - 8 

WORD 

A8 

Enter  the  words  PORE  PRE 

2 

1 - 10 

BETA 

F10.0 

Time  integration  parameter,  the 
recommended  value  is  BETA=0.5 

i)  Undrained  analysis  option 

This  option  is  used  uhen  an  undrained  analysis  is  being  performed 
without  the  use  of  elements  possessing  pore  pressure  degrees  of  freedom. 
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Undrained  analysis 

Card  Columns 

Variable 

Format 

Explanation 

1 1-8 

WORD 

A8 

Enter  the  uord  UNDRAINE 

j)  Axisymmetric  analysis  option 

This  option  is  used  to  obtain  an  ax i symme tr i c analysis.  The  beam 
element  formulation  is  not  affected  by  use  of  this  option.  Integration 
is  performed  over  one  radian;  loads  should  be  calculated  accordingly. 


Axisymmetric  analysis 

Card  Columns  Variable 

Format 

Explanation 

1 1-8  WORD 

A8 

Enter  the  uord  AXISYMME 

k)  Plane  strain  analysis  option 

This  option  is  used  to  obtain  a plane  strain  analysis.  The  beam 
element  formulation  is  not  affected  by  use  of  this  option,  since  it 
assumes  a plane  stress  analysis. 


Pi  ane 

strain  analysis 

Card 

Columns  Variable 

Format 

Exp  1 anat i on 

1 

1 - 8 WORD 

A8 

Enter  the  uords  PL  STRAI 
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1)  Type  of  materials  group 


This  group  specifies  the  material  types  used  in  the  analysis.  For 
details  of  the  available  material  types,  users  should  consult  the 
material  library  section. 


Type 

of  materials 

Card 

Columns  Variable 

Format 

Expl anat i on 

1 

1 - 8 WORD 

A8 

Enter  the  words  TYPE  MAT 

2 

1 - 5 N 

15 

Number  of  different  material  types 

6 - NTYPE(I) 

• 

515 

List  of  material  type  numbers  used 

m)  Type 

of  elements  group 

Th  i s 

group  specifies 

the  element  types  used  in  the  analysis.  For 

details 

of  the  available  element  types  users  should  consult  the  element 

1 i brary 

sect i on . 

Type 

of  elements 

Card 

Columns  Variable 

Format 

Explanation 

1 

1 - 8 WORD 

A8 

Enter  the  words  TYPE  ELE 

2 

1 - 5 N 

15 

Number  of  different  element  types 

6 - NTYPECI ) 

815 

List  of  element  type  numbers  used 

n)  End  of  control  block  card 


This  card  is  required  to  define  the  end  of  the  control  block  input. 
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End  of 

control 

b 1 ock 

Card 

Co  1 umns 

Variable 

Format 

Exp  1 anat i on 

1 

1 - 8 

WORD 

A8 

Enter  the  uords  END  CONT 

A. 6. 2 Mesh  Block 

This  block  contains  the  detailed  information  describing  the  mesh 
geometry,  its  materials  and  boundary  conditions.  This  entire  block 
should  be  omitted  if  a continuation  run  is  being  made  using  the  RESTART 
option. 

a)  Nodal  coordinate  group 

The  coordinates  are  input  for  each  node.  They  may  be  placed  in  any 
order . 


Nodes 


Card 

Columns 

Var i abl e 

Format 

1 

1 - 8 

WORD 

A8 

2 

1 - 5 

N 

15 

6-10 

NCN 

15 

3 

1 - 5 

NNOD 

15 

6 - 

X ( I , NNOD ) 

6F10.0 

Exp  1 anat i on 

Enter  the  word  NODES 
Number  of  nodes  in  this  group 
Number  of  coordinates  per  node 
Node  number 

Nodal  coordinates  1=1, NCN 


3rd  card  repeated  for  a total  of  N times 


b)  Element  connectivity  group 
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The  connectivity  is  input  for  each  element.  The  elements  may  be 
placed  in  any  order.  The  numbering  sequence  for  a given  element  may  be 
found  in  the  element  library  section.  This  group  will  have  to  be  used 
more  than  once  if  more  than  one  element  type  or  more  than  one  material 
is  used. 


El ements 


Card 

Col umns 

Vari abl e 

F ormat 

1 

1 - 

8 

WORD 

A8 

2 

1 - 

5 

N 

15 

6 - 

10 

NETY 

15 

11  - 

15 

NMAT 

I 5 

3 

1 - 

5 

NELE 

15 

6 - 

LC  C I , NELE ) 

2015 

Explanation 

Enter  the  word  ELEMENT 

Number  of  elements  in  this  group 

Element  type  number 

Material  number 

Element  number 

Element  connectivity  1=1, NNE 
NNE  = number  of  nodes  in  element 


3rd  card  repeated  for  a total  of  N times 


c)  Material  properties  group 

The  material  properties  are  input  for  each  material.  The  materials 
may  be  placed  in  any  order.  The  type  of  parameter  to  be  specified  for 
each  material  property  is  shown  in  the  material  library  section.  If 
more  than  one  material  type  is  used,  it  will  be  necessary  to  use  this 
group  more  than  once. 
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Material  properties 


Card 

Co  1 umns 

Variable 

Format 

1 

1 - 

8 

WORD 

A8 

2 

1 - 

5 

N 

I 5 

6 - 

10 

NMTY 

15 

3 

1 - 

5 

NMAT 

I 5 

6 - 

75 

PR ( I , NMAT ) 

7F10.0 

4 

1 - 

PR  Cl , NMAT ) 

8F  10. 0 

Explanation 

Enter  the  word  MATERIAL 
Number  of  materials  in  this  group 
Material  type  number 
Material  number 
Material  properties  for  1=1,7 
Material  properties  for  I=8,NPRM 
NPRM  = number  of  properties  for 
this  material  type 


3rd  and  4th  cards  are  repeated  in  pairs  for  a total  of  N times 


d)  Boundary  conditions  group 

The  boundary  conditions  are  input  in  this  group.  The  total  number 
of  constraints  should  equal  the  value  of  NBC  specified  in  the  control 
block. 


Boundary  conditions 


Card 

Col umns 

Variable 

Format 

1 

1 - 8 

WORD 

A8 

2 

1 - 5 

NSETS 

15 

3 

1 - 5 

N 

15 

6-10 

NDOF 

15 

11  - 20 

ADIS 

F10.0 

4 

1 - 

NNBC(I) 

1615 

Exp  1 anat i on 

Enter  the  word  BOUNDARY 
Number  of  sets  of  constraints 
Number  of  nodes  in  this  set 
Degree  of  freedom  number  being 
constrained  in  this  set 
Applied  displacement  for  this  set 
Node  number  1=1, N 


3rd  and  4th  cards  are  repeated  in  pairs  for  a total  of  NSETS  times 


e)  Mesh  generation  group 
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This  group  may  be  used  to  generate  both  nodal  coordinates  and 


element  connectivity  for  four-  and  eight-noded  quadr i 1 atera 1 elements. 
The  mesh  so  generated  may  be  rectangular  or  radial.  The  first  node  has 
the  smallest  values  of  x and  y as  coordinates.  The  node  numbers  then 
increase  first  along  a line  of  constant  y.  The  elements  are  similarly 
numbered.  This  group  may  be  used  to  generate  all  or  part  of  the  mesh, 
and  may  be  used  more  than  once. 


Mesh  generator 


Card 

Col umns 

Var i abl e 

Format 

Explanation 

1 

1 

- 

8 

WORD 

A8 

Enter  the  uords  MESH  GEN 

2 

1 

5 

MGTY 

15 

Mesh  generator  type  number 
MGTY=1  for  a rectangular  mesh 
MGTY =2  for  a radial  mesh 

6 

- 

10 

NETY 

15 

Element  type  number 

1 1 

- 

15 

NMAT 

15 

Material  number 

3 

1 

- 

5 

NN  1 

15 

First  node  number  generated 

6 

- 

10 

NE  1 

15 

First  element  number  generated 

1 1 

15 

NEX 

I 5 

Number  of  elements  in  the  x or  r 
direction 

16 

— 

20 

NEY 

15 

Number  of  elements  in  the  y or  9 
d i recti  on 

4 

1 

CX(I) 

8F10.0 

Corner  nodal  coordinates  in  the  x 
or  r direction  for  I=1,NEX+1 

5 

1 

CY  ( I ) 

8 F 1 0 . 0 

Corner  nodal  coordinates  in  the  y 
or  0 direction  for  I=1,NEY+1 

f)  End  of  mesh  block  card 


This  card  is  required  to  define  the  end  of  the  mesh  block  input. 
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End  of  mesh  block 


Card  Columns  Variable  Format  Explanation 

1 1-8  WORD  A8  Enter  the  words  END  MESH 


A . 6 . 3 Initial  Stress  Block 

This  block  sets  up  the  initial  stress  conditions  in  the  soil.  If 
zero  stresses  are  desired  then  only  the  end  of  initial  stress  block  card 
is  needed.  The  entire  block  should  be  omitted  if  a continuation  run  is 
being  made  using  the  RESTART  option. 

a)  Initial  stress  variable  group 

Element  stress  and  strain  variables  may  be  explicitly  set  by  using 
this  group.  The  order  of  the  stress  input  depends  on  the  material  type. 
This  order  is  the  same  as  the  stress  output  order  and  may  be  found  in 
the  material  library  section. 
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Initial  stress 

es 

Card 

Co  1 umns 

Variable  Format 

Explanation 

1 

1 - 8 

WORD 

A8 

Enter  the  word  STRESS 

2 

1 - 5 

NSETS 

I 5 

Number  of  initial  condition  sets 

3 

1 - 5 

NE  1 

15 

First  element  number 

6-10 

NE2 

15 

Last  element  number 

11-15 

NI  1 

I 5 

First  integration  point  number 

16  - 20 

N I 2 

I 5 

Last  integration  point  number 

21  - 25 

NSV 

I 5 

Number  of  stress  variables  to  be 
set  at  each  integration  point 

4 

1 - 

SI  G ( I ) 

8F10.0 

Initial  stress  values  for  1=1, NSV 

3rd  and  4th  cards  are 

repeated  in  pairs  for  a total  of  NSETS  times 

o)  Initial  pore  pressure 

group 

Th  i s 

group  is 

used  to 

set  the 

pore  pressure  at  nodal  points.  Pore 

pressures 

may  be 

set  at 

the  integration  points  by  using  the  previous 

group . 

Initial  pore  pressures 

Card 

Col umns 

Variable  Format 

Explanation 

1 

1 - 8 

WORD 

AS 

Enter  the  words  PORE  PRE 

2 

1 - 5 

N 

15 

Number  of  nodes  in  this  group 

3 

1 - 5 

NNOD 

15 

Node  number 

6-15 

U 

F 1 0 . 0 

Pore  pressure 

3rd  card  repeated  for 

a total  of 

N times 

c)  Initial  time  group 


The  total  and  incremental  times  may  be  set  in  this  group. 
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Initial  time 


Card  Columns  Variable  Format  Explanation 

1 1-8  WORD  A8  Enter  the  word  TIME 

2 1-5  T I ME  I F10.0  Incremental  time 

6-10  TIMET  F10.0  Total  time 


d)  Gravity  group 

This  group  is  used  to  simulate  gravity  turn-on  for  material  type  2, 
Cam  clay.  This  material  has  no  strength  at  zero  initial  stress.  Hence, 
it  is  not  possible  to  use  the  normal  procedure  of  gravity  turn  on.  The 
stresses  are  calculated  at  each  integration  point  as  shown  below,  where 
the  depth,  d,  is  taken  to  be  the  negative  of  the  y coordinate  at  that 
point.  The  pore  pressure,  u,  is  given  by 


u = drw  + qw  (A. 20) 

The  effective  stresses  are  given  by 

ay/=dyt  + q-t-u  (A. 21) 
ax/=koay'  (A. 22) 
Oz'  = ko^y'  (A. 23) 

where  the  constants  qw  and  q t may  be  used  to  simulate  a linear  change  of 
stress  with  depth  for  any  water  table  and  overburden  pressure.  The 
preconsolidation  pressure  is  given  by 


Pc' 


OCR 


q2  1 

p'  + 

p'm2J 


(A. 24) 
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Gravity 

Card  Columns 

Variable 

Format 

Exp  1 anat i on 

1 1-8 

WORD 

A8 

Enter  the  word  GRAVITY 

2 1-5 

NSETS 

I 5 

Number  of  sets  of  elements  for 
which  gravity  properties  are  input 

3 1-5 

N E 1 

I 5 

First  element  number 

6-10 

NE2 

I 5 

Last  element  number 

11  - 20 

PR  G ( 1 ) 

F10.0 

y t unit  we  i glit  of  soil 

21  - 30 

PRGC2) 

F10.0 

■y  w unit  weight  of  water 

31  - 40 

PR G ( 3) 

F 1 0 . 0 

k0  the  coefficient  of  earth 
pressure  at  rest 

O 

LO 

1 

TT 

PRG ( 4) 

F10.0 

OCR  the  overconso 1 i dat i on  ratio 

51  - 60 

PRG( 5) 

F 10. 0 

qt  the  overburden  constant 

61  - 70 

PRG ( 6 ) 

F10.0 

qw  the  water  table  constant 

3rd  card  repeated  for  a 

total  of 

NSETS  times 

e)  End  of  initial 

stress  block  card 

This  card  is 

required 

to  define 

the  end  of  the  initial  stress  input. 

End  of  initial 

stress  block 

Card  Columns 

Variable 

Format 

Explanation 

1 1-8 

WORD 

A8 

Enter  the  words  END  INIT 

A.  6. 4 Load  Block 

This  block  is  used  to  apply  loads.  The  entire  block  may  be  repeated 
so  that  incremental  load  problems  may  be  solved.  The  maximum  number  of 
load  increments  applied  in  one  run  is  controlled  by  the  value  of  NLLI, 
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which  is  input  in  the  control  block. 


Execution  will  stop  before  NLL  I 


increments  if  insufficient  load  increment  data  is  provided. 

The  incremental  load  vector  is  initially  set  to  zero,  but  it  is  not 
reset  to  zero  between  increments.  Thus,  nodal  loads  applied  in  one 
increment  will  be  applied  again  in  the  next  increment  unless  PINC  is  set 
to  zero  in  the  proportional  increment  group. 

a)  Nodal  load  group 

This  group  may  be  used  to  apply  nodal  loads.  These  incremental 
loads  are  added  to  the  existing  values  in  the  incremental  load  vector. 


Nodal 

1 oads 

Card 

Co  1 umns 

Variable 

Format 

Exp  1 anat i on 

1 

1 - 8 

WORD 

A8 

Enter  the  words  NODAL  LO 

2 

1 - 5 

N 

15 

Number  of  nodes  in  this  group 

6-10 

NDFN 

I 5 

Number  of  degrees  of  freedom 
at  this  node 

3 

1 - 5 

NNOD 

15 

Node  number 

6 - 

A ( I ) 

7F10.0 

Nodal  loads  1=1, NDFN 

3rd  card  repeated  for  a 

total  of 

N times 

b)  Body  force  group 

This  group  is  used  to  input  body  forces.  These  forces  are  converted 
into  equivalent  nodal  loads  which  are  then  added  to  the  existing 
incremental  load  vector. 
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Body  forces 

Card  Columns 

Variable 

Format 

Exp  1 anat i on 

1 1-8 

WORD 

A8 

Enter  the  words  BODY  FOR 

2 1-5 

NSETS 

15 

Number  of  sets  of  elements  for 
which  body  forces  are  input 

3 1-5 

NE  1 

15 

First  element  number 

6-10 

NE2 

I 5 

Last  element  number 

11  - 20 

FX 

F10.0 

x force  per  unit  volume 

21  - 30 

FY 

F10.0 

y force  per  unit  volume 

31  - ‘10 

FZ 

F10.0 

z force  per  unit  volume 

3rd  card  repeated  for  a total  of 

NSETS  times 

c)  Proportional 

increment  group 

This  group 

may  be  used 

for  two  purposes;  to  scale  the  existing 

incremental  load 

vector,  and 

secondly,  to  apply  repeated  load  increments 

without  repeating  the  entire  load 

block.  The  incremental  load  vector 

may  be  reset  to 

zero  by  using  this 

group  with  PINC^O.O. 

Proport i onal 

i ncrement 

Card  Columns 

Variable 

F ormat 

Exp  1 anat i on 

1 1-8 

WORD 

A8 

Enter  the  word  PROPORTI 

2 1-5 

NPINC 

15 

The  number  of  load  increments  for 
which  PINC  is  to  be  applied 

6-15 

PINC 

F10.0 

The  incremental  load  vector  is 
scaled  by  this  factor 

d>  Time  increment  group 
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This  group  is  used  to  control  the  time  increment  for  problems  which 


have  pore  pressure  degrees  of  freedom.  After  each  increment  the  total 
time  is  incremented  by  adding  the  incremental  time;  however,  the 
incremental  time  is  not  reset  to  zero.  Thus,  if  the  incremental  time  is 
set  in  the  initial  stress  block  and  this  group  is  not  used  the  time 
increment  will  remain  constant.  For  most  consolidation  problems  this  is 
undesirable,  since  simulating  conditions  from  undrained  to  drained 
usually  requires  about  a 1000  to  10000  fold  increase  in  time.  The 
following  two  schemes  are  provided  for  altering  this  time  step 


1.  The  manual  method  is  implemented  by  explicitly  setting  the  time 
increment  for  each  step.  This  is  done  by  setting  NTINC=0  and 
using  this  group  in  each  repeated  load  block. 


2.  The  automatic  method  is  to  specify  how  many  increments  are 
desired  to  achieve  a tenfold  increase  in  time.  For  example, 
TINC=5.0  allows  the  time  to  increase  by  a factor  of  ten  every 
five  increments.  The  time  increments  will  appear  evenly  spaced 
when  plotted  on  a logarithmic  scale.  The  incremental  time  is 

computed  using  the  following  formula: 


ti 


' 1/TINC 
t0  10 


1 


(A. 25) 


where  to  is  the  total  time  at  the  start  of  the  current  increment. 
Note  that  the  value  of  t0  must  not  be  zero. 
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Time  increment 

Card 

Co  1 umns 

Variable 

Format 

Exp  1 anat i on 

1 

1 - 8 

WORD 

A8 

Enter  the  words  TIME  INC 

2 

1 - 5 

NTINC 

15 

Number  of  time  increments 
NTINC=0  for  absolute  increment 

6-15 

TINC 

F10.0 

Number  of  increments  for  a tenfold 
increment  in  total  time 
If  NTINC=0  absolute  increment 

e)  Change 

material 

group 

Th  i s 

group  is 

used  to  change 

the  material  properties  of  an  existing 

material . 

With  careful  choice  of 

properties  it  may  be  used  to  simulate 

the  construction  of  embankments  or 

an  excavation. 

Change 

mater i al 

Card 

Col umns 

Variable 

Format 

Exp  1 anat i on 

1 

1 - 8 

WORD 

A8 

Enter  the  words  CHANGE  M 

2 

1 - 5 

NMAT 

15 

Material  number 

6-10 

NPRM 

15 

Number  of  properties  for  this 
material 

3 

1 - 

PR (I, NMAT) 

8F10. 

0 New  properties  1=1, NPRM 

f)  End  of  load  block  card 

This  card  is  required  to  define  the  end  of  the  load  block.  This 
block  should  be  repeated  if  more  than  one  load  increment  is  desired, 
unless  NPINC  or  NTINC  is  greater  than  one. 
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End  of  load  block 


Card  Columns  Variable  Format  Explanation 

1 1-8  WORD  AS  Enter  the  words  END  LOAD 


A . 7 SAMPLE  INPUT 

The  input  for  the  simple  problem  of  one-dimensional  consolidation  of 
an  elastic  layer  is  given  below.  A mesh  consisting  of  five  L3P2 
elements  is  used.  The  displacement  at  node  eleven  is  fixed,  and  the 
pore  pressure  at  node  one  is  fixed.  A uniform  load  is  applied  at  node 
one  during  the  first  increment.  A total  of  fourteen  time  increments  are 
used . 


TITLE 

ONE  DIMENSIONAL  CONSOLIDATION  OF  AN  ELASTIC  LAYER 
SIZE 

11  5 1 2 

PORE  PRESSURE 
0.5 

TYPE  ELEMENT 
1 7 

TYPE  MATERIAL 
1 4 

INCREMENTS 

1 14 
END  CONTROL 
NODES 

1 1 1 

1 0.0 

2 -10.0 

3 -20.0 

4 -30.0 

5 -40.0 

6 -50.0 

7 -60.0 

8 -70.0 

9 -80.0 
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10 

o 

o 

1 

1 1 

-100.0 

ELEMENTS 

5 

7 1 

1 

3 1 

2 

2 

5 3 

4 

3 

7 5 

6 

4 

9 7 

8 

5 

1 1 9 

10 

MATERIALS 

1 

4 

1 

0.0 

0.0 

o 

o 

62.4 

0.0433 

14-1000 

.0 

0.0 

BOUNDARY  CONDITIONS 
2 

1 1 

1 1 

1 2 

1 

END  MESH 
END  INITIAL 
NODAL  LOADS 
1 1 

1 -1.0 

TIME  INCREMENT 

0 0.63 
END  LOAD 

PROPORTIONAL  INCREMENT 

1 0.0 

TIME  INCREMENT 

13  0.2 

END  LOAD 
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Appendix  B 


SUBROUTINE  LISTINGS 

This  appendix  includes  listings  of  FORTRAN  subroutines  similar  to 
those  used  in  PEPCO.  They  are 

1.  DMAT , a subroutine  for  calculating  the  stress-strain  matrix  at  an 
integration  point  for  problems  of  axial  symmetry  or  plane  strain. 

2.  UPDATE,  a subroutine  for  updating  the  Cam  clay  yield  surface  at 
an  integration  point  after  the  incremental  stresses  have  been 
recovered . 

3.  FLOW,  a subroutine  for  expanding  a Q8  element  to  a Q8P4  element 
in  order  to  include  the  effect  of  consolidation. 

4.  LARGE,  a subroutine  for  calculating  the  initial  stress  stiffness 
matrix  which  is  to  be  added  to  the  element  stiffness  matrix  in 
order  to  include  the  effect  of  large  strains. 

5.  ROTATE,  a subroutine  to  calculate  the  true  stress  increment  as  a 
result  of  rotations. 
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SUBROUTINE  DMAT ( D , SVT , PR , DEPTH ) 

C 

C SUBROUTINE  CALCULATES  THE  STRESS-STRAIN  MATRIX  FOR  CAM  CLAY 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  SVTC  1 ) ,PRC 1 ) , DC4, 4) ,A1 C4) ,A2C4) 

PC=SVT( 10) 

P-SVT ( 1 1 ) 

Q-SVT ( 12) 

E -SVT (13) 

C 

C CALCULATE  THE  ELASTIC  STRESS-STRAIN  MATRIX 
C 

BMOD=( 1 .+E)*P/PR( 1 ) 

SU=PR(6)+PR(7)*DEPTH 
SMOD=PR( 5)*SU 
Dl=(3.#BH0D+4.#SM0D)/3. 

D2=C3.*BM0D-2.#SM0D)/3. 

DC  1 , 1)=D1 
DC  1 ,2)=D2 
D C 1 , 3) =D2 
D C 1 > 4 ) =0 . 

D ( 2 , 1 )=D2 
D C 2 > 2 ) =D  1 
D ( 2 , 3)=D2 
D(2,4)=0. 

DC3, 1 ) =02 
DC3,2)=D2 
D(3, 3) =D  1 
D C 3, 4 ) =0  . 

D ( 4 , 1 ) =0 . 

D(4,2)=0. 

D ( 4, 3) =0  . 

D ( 4 , 4 ) -SMO  D 
C 

C CHECK  TO  SEE  IF  THIS  POINT  LIES  ON  THE  YIELD  SURFACE 
C 

AM-PR  C 4 ) 

PO=P+Q*Q/(P*AM*AM) 

IF  (PO. LT. 0. 999*PC)  RETURN 
C 

C CALCULATE  THE  ELASTO-PLAST I C STRESS-STRAIN  MATRIX 
C 

Dl=(2.*P-PC)/3. 

D2=3./(AM*AM) 

A1C1)=D1+CSVT(5)-P)*D2 
A1  (2)=D1  + (SVTC6)-P)*D2 
A1 C3)=D1+(SVT(7)-P)*D2 
A1C4)=2.*SVTC8)*D2 
DO  22  1=1,3 
TEMP=0 . 

DO  21  J= 1 , 3 

21  TEMP=TEMP+DCI , J ) *A  1 (J) 

22  A2 C I ) =TEMP 
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A2(4)=SMOD*A1  (4) 

TEMP=0 . 

DO  23  1=1,4 

23  T EMP  = T EMP+A 1 (I )*A2(I ) 

TEMP=TEMP+P*PC#(1 .+E)#(2.*P-PC)/(PR(2)-PR( 1)) 
DO  24  1=1,4 

DO  24  J=1 , 4 

24  D< I , J)=D( I , J)-A2(I )#A2( J)/TEMP 
RETURN 

END 


SUBROUTINE  UPDATE (SV I , SVT , PR ) 

C 

C SUBROUTINE  TO  UPDATE  CAM  CLAY  DEPENDENT  STRESS  QUANTITIES 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  SVI(1),SVT(1),PR(1) 

AK  = PR ( 1 ) 

AL  = PR ( 2 ) 

EO=PR ( 3) 

AM=PR (4) 

SX=SVT ( 5) 

SY=SVT (6) 

SZ=SVT ( 7 ) 

TX Y = SVT ( 8 ) 

PC  = SVT(  10) 

P=(SX+SY+SZ)/3. 

Q=DSQRT(SX*(SX-SY)+SY*(SY-SZ)+SZ*(SZ-SX)+3.*TXY*TXY) 

IF  (P.LT.0.01)  P=0.01 

PO=P+Q*Q/(P*AM*AM) 

IF  (PO.GT.PC)  PC=PO 

E=EO-AK*DLOG(P)-(AL-AK)*D LOG (PC/2. ) 

SVT(10)=PC 
SVT ( 1 1 ) =P 
SVT ( 1 2 ) =Q 
SV  T ( 1 3) =E 
RETURN 
END 
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SUBROUTINE  FLOW(SE,AE,PP,PR,F4,B4,B8,DETJ4,DETJ8,WGHT,TIMEI, 
2BETA , N , N2 ) 

C 

C SUBROUTINE  TO  EXPAND  THE  ELEMENT  STIFFNESS  MATRIX  TO  INCLUDE  A PORE 
C PRESSURE  DEGREE  OF  FREEDOM  AT  THE  CORNER  NODES,  AND  ALSO  TO 
C CALCULATE  THE  CORRESPONDING  RIGHT  HAND  SIDE  VECTOR 
C 

IMPLICIT  REAL#8(A-H,0-Z) 

DIMENSION  SE(1),AE(1),PP(N),PR(7),F4(N),B8(4,N2),E(3,8),PE(3,8), 
2ETPE(8,8),RN(16),B4(4,8) 

GAMMAW=PR(4) 

PERMX-PR ( 5 ) /GAMMAU 
PERMY=PR(6 ) /GAMMAN 
DO  1 1=1, N 
E ( 1 , I ) =B4 ( 1,1) 

E ( 2 , I ) =B4 ( 2 , I+N ) 

1 E ( 3, I )=0. 0 
DO  2 1 = 1 ,N 

PE(1,I)=PERMX*E(1,I) 

PE(2,I)=PERMY*E(2»I ) 

2 PE ( 3, I )=PERMX*E(3, I ) 

DO  3 1=1, N 

DO  3 J=1 ,N 
T EMP  = 0 . 0 
DO  4 K= 1 , 3 

4 TEMP=TEMP+E(K, I )*PE(K, J) 

3 ETPE ( I , J) =TEMP*T IMEI *DETJ4*UGHT 
DO  5 1 = 1 , N 2 

5 RN( I )=B8( 1 , I )+B8(2, I )+B8(3, I ) 

N2T  = (N2+ 1 ) *N2/2 

L=N2T 

DO  6 J= 1 , N 
L=L+J- 1 
DO  6 1=1, N2 
L = L+ 1 

6 SE(L)=SE(L)+RN(I ) *F 4 ( J ) *DET J8*UGHT 
L = N2T 

DO  7 J=1 ,N 
L=L+N2 
DO  7 1=1, J 
L =L+  1 

7 SE(L)=SE(L)-ETPE(I,J)*BETA 
DO  8 1=1, N 

K=I+N2 
TEMP=0 . 0 
DO  9 J=  1 , N 

9 TEMP=TEMP+ETPE C I , J)#PP(J) 

G I N=0 . 

DO  10  J= 1 , 3 

10  GIN=GIN+PR(J)*PE(J, I ) 

GIN=GIN*TIMEI*DETJ4*WGHT*PR(4) 

8 AE(K)=AE(K) +T EMP+G I N 
RETURN 

END 
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SUBROUTINE  LARGE ( SE , S I G , PFS , PFT , DET J , UG HT , N ) 

C 

C SUBROUTINE  TO  CALCULATE  THE  INITIAL  STIFFNESS  MATRIX  FOR  A LARGE 
C DEFORMATION  ANALYSIS  FOR  A N NODED  ISOPARAMETRIC  QUADRILATERAL 
C ELEMENT  IN  PLANE  STRAIN 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  SE(1),SIG(4),PFS(N),PFT(N),SSS(8,8),STT(8,8),SST(8,8) 

XX J=-S I G ( 1 )*DETJ*UGHT 

YYJ=~SIG(2)*DETJ*  LIGHT 

XYJ=-SIGH)*DETJ*I4GHT 

HD  I F -0 . 5* ( YYJ-XXJ  ) 

H S U M = 0 . 5*(XXJ+YYJ) 

DO  1 1=1, N 
DO  1 J=1,N 

SSS(I ,J)=PFS(I )*PFS(J) 

ST  T ( I , J)=PFT(I )*PFT(J) 

1 SST  C I , J)=PFT(I )*PFS(J) 

NT  = ( N+ 1 )*N/2 

L = 0 

DO  2 J = 1 > N 
DO  2 1 = 1 , J 
L = L+  1 

SE(L)=SE(L)+HDIF*STT(I , J ) -XX J *SSS ( I , J) 

K = L + NT  + J*N 

2 SE (K)=SE (K)-HDI F*SSS( I , J)-YYJ*STT (I , J) 

L = NT 

DO  3 J=1 ,N 
L = L-+ J-  1 
DO  3 1=1, N 
L = L+  1 

3 SE(L)=SE(L)-HSUM^SST( I , J ) -XY J* ( SSS ( I , J )+STT ( I , J ) ) 

RETURN 

END 


SUBROUTINE  ROTATE (SI GT,SI GI , EPSI , B, DX,N,N2) 

C 

C SUBROUTINE  TO  MAKE  A CORRECTION  TO  THE  INCREMENTAL  STRESSES  DUE  TO 
C THE  EFFECT  OF  ROTATIONS,  FOR  A N NODED  ISOPARAMETRIC  QUADRILATERAL 
C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  SIGT(4),SIGI(4),EPSI(4),B(4,N2),DX(N2) 

PU Y=0 . 0 
PVX=0. 0 
DO  1 1=1, N 

PUY=PUY+B(4, I )*DX(I ) 

J = I + N 

1 PVX=PVX+B(4,J)*DX(J) 

R=0 . 5*( PVX-PUY ) 

E=EPSI ( 1 )+EPSI (2)+EPSI (3) 
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SIGI ( 1 )=SIGI ( 1 )+SIGT( 1 )*E+2.*R*SIGT(4) 

SIGI (2)=SIGI (2)+SIGT(2)*E~2.*R*SIGT(4) 
SIGI(3)=SIGI(3)+SIGT(3)*E 

SIGI (4)=SIGI (4)+SIGT(4)*E~R*(SIGT( 1 )-SIGT(2)  ) 

RETURN 

END 
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Appendix  C 


TUNNEL  INPUT  DECK  FOR  PEPCO 


This  appendix  provides  a listing 
tunnel  consolidation  analysis.  The 
Chapter  V I . 


of  the  input  deck  for  a typical 
finite  element  mesh  is  given  in 


TITLE 

TUNNEL  IN  BOSTON  BLUE  CLAY,  DIAMETER  = 8M,  DEPTH  TO  CROHN  = 9.5M 
TITLE 

LINER  ACTIVATED  AFTER  3.6  INCHES  OF  DISPLACEMENT  AT  THE  CROWN 
SIZE 

263  74  2 88 

PL  STRAIN 
PORE  PRESSURE 
0.5 

TYPE  ELEMENT 

2 3 4 

TYPE  MATERIAL 


2 

1 2 

INCREMENTS 

1 

63 

END  CONTROL 

NODES 

263 

2 

1 

0.0000 

-32.8000 

2 

2.2396 

-33.0207 

3 

4.3932 

-33. 6738 

4 

6.3780 

-34.7349 

5 

8. 1177 

-36. 1623 

6 

9.5451 

-37.9020 

7 

10.6062 

-39.8867 

8 

1 1 .2593 

-42.0404 

9 

1 1 . 4800 

-44.2800 

10 

1 1 .2593 

-46.5196 

1 1 

10.6062 

-48.6732 

12 

9.5451 

-50.6579 

13 

8.1177 

-52.3977 

14 

6.3780 

-53.8251 
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15 

4.3932 

-54.8862 

16 

2.2396 

-55.5392 

17 

0.0000 

-55.7600 

18 

0.0000 

-31 . 9800 

19 

4.7071 

-32.9161 

20 

8.6976 

-35.5827 

21 

1 1 . 3636 

-39.5729 

22 

12.3000 

-44.2800 

23 

1 1 . 3636 

-48.9871 

24 

8.6976 

-52.9772 

25 

4.7071 

-55.6435 

26 

0.0000 

-56.5800 

27 

0.0000 

-31 . 1600 

28 

2.5597 

-31 .4122 

29 

5.0207 

-32. 1588 

30 

7.2891 

-33. 3710 

31 

9.2771 

-35.0028 

32 

10.9090 

-36.9908 

33 

12.1212 

-39.2593 

34 

12.8678 

-41  . 7203 

35 

13. 1200 

-44.2800 

36 

12.8678 

-46.8397 

37 

12.1212 

-49.3007 

38 

10.9090 

-51  . 5691 

39 

9.2771 

-53.5571 

40 

7.2891 

-55. 1889 

41 

5.0207 

-56.4012 

42 

2.5597 

-57. 1477 

43 

0.0000 

-57.4000 

44 

0.0000 

-29.5200 

45 

5.6485 

-30.6434 

46 

10.4370 

-33.8430 

47 

13.6366 

-38.6315 

48 

14.7600 

-44.2800 

49 

13.6366 

-49.9285 

50 

10.4370 

-54.7170 

51 

5.6485 

-57.9162 

52 

0.0000 

-59.0400 

53 

0.0000 

-27.8800 

54 

3.  1996 

-28. 1952 

55 

6.2760 

-29. 1284 

56 

9.1115 

-30.6440 

57 

1 1 . 5964 

-32.6835 

58 

13.6359 

-35. 1688 

59 

15. 1516 

-38.0040 

60 

16.0848 

-41  . 0804 

61 

16.4000 

-44.2800 

62 

16.0848 

-47.4796 

63 

15.1516 

-50.5559 

64 

13. 6359 

-53. 3912 

65 

1 1 . 5964 

-55.8764 

66 

9.1115 

-57.9159 

67 

6.2760 

-59.4316 

68 

3. 1996 

-60.3648 
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69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

1 1 1 

112 

113 

114 

115 

116 

117 

118 

119 

120 

121 

122 


0.0000 

-60.6800 

0.0000 

-25.4200 

7.2173 

-26.8557 

13. 3362 

-30.9438 

17.4243 

-37.0627 

18.8600 

-44.2800 

17.4243 

-51 . 4973 

13.3362 

-57.6161 

7.2173 

-61 . 7043 

0.0000 

-63. 1400 

0.0000 

-22.9600 

4.1594 

-23. 3697 

8.  1587 

-24.5829 

1 1.8447 

-26.5532 

15.0755 

-29.2045 

17.7271 

-32.4353 

19.6971 

-36. 1213 

20.9103 

-40.  1 206 

21 . 3200 

-44.2800 

20.9103 

-48.4394 

19.6971 

-52 .4387 

17.7271 

-56. 1247 

15.0755 

-59.3555 

1 1.8447 

-62.0067 

8. 1587 

-63. 9771 

4.1594 

-65.  1903 

0.0000 

-65.6000 

0.0000 

-19.6800 

9.4139 

-21 . 5525 

17.3948 

-26.8852 

22.7274 

-34.8661 

24.6000 

-44.2800 

22.7274 

-53.6939 

17.3948 

-61 . 6743 

9.4139 

-67.0074 

0.0000 

-68.8800 

0.0000 

-16.4000 

5.4392 

-16.9356 

10.6692 

-18.5222 

15.4895 

-2  1 . 0986 

19.7141 

-24.5659 

23.1814 

-28.7908 

25.7578 

-33. 6108 

27. 3444 

-38.8408 

27.8800 

-44.2800 

27.3444 

-49.7192 

25.7578 

-54.9492 

23. 1814 

-59.7691 

19.7141 

-63.9941 

15.4895 

-67.4614 

10.6692 

-70.0378 

5.4392 

-71 . 6244 

0.0000 

-72. 1600 

0.0000 

-12.3000 
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123 

12.2383 

-14.7344 

124 

22.6133 

-21.6667 

125 

29.5456 

-32.0417 

126 

31 . 9800 

-44.2800 

127 

29. 5456 

-56.5183 

128 

22.6133 

-66.8933 

129 

12.2383 

-73.8256 

130 

0.0000 

-76.2600 

131 

0.0000 

-8.2000 

132 

7.0389 

-8.8934 

133 

13.8072 

-10.9463 

134 

20.0451 

-14.2805 

135 

25.5125 

-18.7675 

136 

29.9995 

-24.2349 

137 

33. 3336 

-30.4728 

138 

35.3866 

-37.241 1 

139 

36.0800 

-44.2800 

140 

35. 3866 

-51 . 3189 

141 

33.3336 

-58.0872 

142 

29.9995 

-64.3250 

143 

25.5125 

-69.7925 

144 

20.0451 

-74.2792 

145 

13.8072 

-77.6137 

146 

7.0389 

-79.6666 

147 

0.0000 

-80.3600 

148 

0.0000 

-4.  1000 

149 

14.2837 

-5.4733 

150 

27.5162 

-9. 3838 

151 

34.8962 

-16.7637 

152 

38.8070 

-29.9966 

153 

40.  1800 

-44.2800 

154 

38.8070 

-58.5637 

155 

34.8962 

-71.7962 

156 

27.5162 

-79.1762 

157 

14.2837 

-83.0866 

158 

0.0000 

-84.4600 

159 

0.0000 

0.0000 

160 

7.3800 

0.0000 

161 

14.7600 

0.0000 

162 

22. 1400 

0.0000 

163 

29.5200 

0.0000 

164 

36.9000 

0.0000 

165 

44.2800 

0.0000 

166 

44.2800 

-7.3800 

167 

44.2800 

-14.7600 

168 

44.2800 

-22. 1400 

169 

44.2800 

-29.5200 

170 

44.2800 

-36.9000 

171 

44.2800 

-44.2800 

172 

44.2800 

-51.6600 

173 

44.2800 

-59.0400 

174 

44.2800 

-66.4200 

175 

44.2800 

-73.8000 

176 

44.2800 

-81 . 1800 

177 

178 

179 

180 

181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

211 

212 

213 

214 

215 

216 

217 

218 

219 

220 

221 

222 

223 

224 

225 

226 

227 

228 

229 

230 


44.2800 

-88.5600 

36.9000 

-88.5600 

29.5200 

-88.5600 

22.  1400 

-88.5600 

14.7600 

-88.5600 

7.3800 

-88.5600 

0.0000 

-88.5600 

51  . 6600 

0.0000 

51  . 6600 

-14.7600 

51 . 6600 

-29. 5200 

51.6600 

-44.2800 

51.6600 

-59.0400 

51 . 6600 

-73. 8000 

51 . 6600 

-88.5600 

59.0400 

0.0000 

59.0400 

-7. 3800 

59.0400 

-14.7600 

59.0400 

-22.1400 

59.0400 

-29.5200 

59.0400 

-36.9000 

59.0400 

-44.2800 

59.0400 

-51 . 6600 

59.0400 

-59.0400 

59.0400 

-66.4200 

59.0400 

-73. 8000 

59.0400 

-81 . 1800 

59.0400 

-88.5600 

68.8800 

0.0000 

68.8800 

-14.7600 

68.8800 

-29.5200 

68.8800 

-44.2800 

68.8S00 

-59.0400 

68.8800 

-73. 8000 

68.8800 

-88.5600 

78.7200 

0.0000 

78.7200 

-7.3800 

78.7200 

-14.7600 

78.7200 

-22.  1 400 

78.7200 

-29.5200 

78.7200 

-36.9000 

78.7200 

-44.2800 

78.7200 

-51 . 6600 

78.7200 

-59.0400 

78.7200 

-66.4200 

78.7200 

-73. 8000 

78.7200 

-81 . 1800 

78.7200 

-88.5600 

91 . 8400 
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Appendix  D 

REPORT  OF  NEW  TECHNOLOGY 

The  work  performed  under  this  contract  has  lead  to  no  new  technolog 
cal  inventions.  Conclusions  and  recommendations  regarding  various  types 
of  equipment  and  procedures,  design  parameters,  and  soi 1 /structure  inter 
action  are  intended  to  expand  and  improve  the  state-of-the-art  of  tunnel 
design  and  construction  on  soft  ground. 
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